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ABSTRACT 


Large  Static  and  Dynamic  Deformation  of  Beams: 
The  Inverse  Problem.  (August  1996) 

Matilda  Wilson  McVay,  B.S.,  Colorado  School  of  Mines; 
M.S.,  Texas  A&M  University 
Chair  of  Advisory  Committee;  Dr.  John  L.  Junkins 


The  primary  objective  of  this  study  is  to  establish  a  methodology  for  validating  a 
dynamical  model  of  a  flexible  two  beam  system  undergoing  large  deformations. 

This  dissertation  presents  a  systematic  approach  for  using  a  series  of  experiments  to 
estimate  mathematical  model  parameters  and  correct  them  to  match  the  measured  response. 
The  study  of  any  dynamical  system  has  a  natural  partition  between  the  kinetic  energy  terms 
and  the  potential  energy  terms.  By  using  this  natural  partition  to  design  a  sequence  of 
experiments,  it  is  shown  that  the  number  of  unknown  parameters  affecting  the  system 
response  for  any  given  experiment  is  greatly  reduced.  First,  a  set  of  static  deformations  are 
used  to  determine  beam  parameters  such  as  the  stiffness  coefficients  and  allow  modeling  of 
nonlinear  effects.  Then,  free  response  experiments  are  used  to  determine  some  motion 
parameters  such  as  the  mass  per  unit  length  of  each  beam  and  the  parameters  associated 
with  natural  environmental  forces  such  as  friction  effects.  This  separation  of  static  and  free 
response  measurements  allows  the  recovery  of  model  parameters  without  being  corrupted 
by  other  forced  system  model  errors  such  as  joint  dynamics  and  motor  modeling  which  are 
present  in  a  full  dynamic  response.  A  set  of  dynamic  forced  response  experiments  are  used 
to  determine  motor  parameters  which  model  the  inputs  to  the  structure.  Appropriate 
statistical  estimation  methods  are  utilized  to  forward  propagate  a  priori  and  measurement 
covariance  estimates  through  the  sequence  of  nonlinear  estimation  processes. 


Before  the  parameters  can  be  updated,  a  novel  mathematical  model  of  the  system  is 
developed  and  verified  based  on  an  arc  length  approach  to  beam  deflections.  This  model 
accounts  for  large  deformation  affects  such  as  foreshortening  and  is  used  in  the  static 
analysis.  A  nonlinear  finite  element  model  is  developed  to  allow  modeling  of  the  free 
vibration  and  forced  vibration  response.  This  model  also  allows  large  deformations  of  the 
beam  system  by  defining  the  deformation  of  a  point  in  terms  of  an  axial  displacement  and  a 
transverse  displacement.  Although  it  is  a  computationally  intensive  program,  it  yields 
accurate  results  and  is  well  suited  to  the  estimation  of  parameters  affecting  the  dynamic 
motion. 
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CHAPTER  I 


INTRODUCTION 


1.1  Motivation 

Most  of  the  parameters  involved  in  a  mathematical  model  of  a  dynamical  system 
can  be  naturally  grouped  into  three  major  subsets.  The  first  logical  subset  would  include 
parameters  affecting  the  storage  of  potential  energy  which  can  be  estimated  from  a  family 
of  static  measurements  on  the  system  displacements.  The  second  subset  would  include 
inertia  type  parameters  which  affect  the  free  response  of  the  system  to  prescribed  initial 
conditions.  These  could  be  found  by  deforming  the  structure  to  a  measured  static  shape 
and  releasing  the  constraining  force  and  measuring  the  free  response  of  the  system.  The 
third  subset  would  contain  the  actuator  model  parameters  which  could  be  found  from  a  set 
of  forced  response  experiments.  This  parameterization  reduces  the  number  of  poorly 
known  parameters  to  be  determined  from  each  identification  experiment  and  typically 
improves  observability.  For  high  dimensioned  nonlinear  systems,  the  approach  would 
seem  likely  to  greatly  enhance  convergence.  These  qualitative  ideas  are  the  basis  for  the 
methodology  studied  in  this  dissertation. 

These  ideas  are  studied  in  the  context  of  large  deformations  of  a  two  beam  system 
that  is  similar  to  the  flexible  manipulator  arm  developed  at  the  United  States  Air  Force 
(USAF)  Phillips  Laboratory  called  the  Planar  Articulating  Controls  Experiment  (PACE)'. 
This  structure,  illustrated  in  Figure  1,  consists  of  two  flexible  beams  connected  to  each 
other  by  an  elbow  joint  and  mounted  to  an  air  bearing  table  by  a  shoulder  hub.  The 
primary  beam  material  properties  and  the  tip  and  elbow  mass  properties  are  listed  in  Table 
1.  The  structure  is  constrained  to  move  horizontally  along  a  polished  granite  table.  Air 
bearings  are  used  to  greatly  reduce  die  friction  forces  between  the  table  and  the  support 


This  dissertation  follows  the  style  and  format  of  the  AIAA  Journal  of  Guidance,  Control 
and  Dynamics. 


2 


•  plates.  PACE  is  designed  to  experimentally  test  control  theories  and  validate  nonlinear 

structural  modeling  techniques.  A  similar  stmcture  is  being  developed  in  the  Dynamics 
and  Controls  Laboratory  of  the  Aerospace  Engineering  Department  at  Texas  A&M 
University. 


Figure  1.  Planar  Articulating  Controls  Experimental  test  article 


1  Table  1.  Material  properties  for  PACE 

Beam  1 

Beam  2 

L=0.776  m 

L=0.714m 

EI=1 1.413  N-m2 

EI=1 1.275  N-m2 

p  =0.532  kg/m 

p  =0.530  kg/m 

Elbow  mass=4.280  kg 

Tip  mass=1.038  kg 

Elbow  inertia=122.982E-4  kg-m^ 

Tip  inertia=15.244E-4  kg-m^ 

1.2  Dissertation  Outline 


This  dissertation  begins  in  Chapter  n  with  a  review  of  the  literature  for  the 
modeling  of  large  deformations,  with  a  specific  focus  on  two  beam  systems.  The  next 
three  chapters  correspond  to  three  divisions  of  the  system  identification  process  into  the 
static,  free  vibration,  and  forced  response  analysis.  Each  chapter  is  divided  into  two  main 
sections.  The  first  section  develops  the  mathematical  model  for  simulating  the  response 
and  the  second  section  develops  the  parameter  estimation  technique.  Some  highlights  of 
each  chapter  are  discussed.  Chapter  HI  develops  a  new  approach  to  the  nonlinear  static 
beam  flexure  problem.  This  model  is  used  in  estimating  the  beam  parameters  that  affect  the 
potential  energy  by  using  a  series  of  simple  static  experiments  designed  to  deforni  the 
system  using  a  wire  constraint  force.  An  unusual  application  of  the  implicit  function 
theorem  is  used  in  this  recovery  process.  Chapter  IV  develops  and  validates  a  nonlinear 
finite  element  model  for  the  dynamic  analysis  of  the  system.  This  model  is  used  to 
estimate  the  parameters  contributing  to  the  free  response  of  the  system  upon  release  of  the 
force  applied  by  the  wire.  Chapter  V  develops  the  additions  to  the  mathematical  model  for 
including  a  motor  in  the  dynamic  response  and  recovers  the  motor  parameters.  A  method 
is  presented  for  propagating  the  estimated  parameters  and  associated  errors  forward 
through  each  level  of  parameter  estimation.  Finally,  Chapter  VI  presents  the  conclusions 
of  the  system  modeling  and  parameter  estimation  technique  studied  in  this  dissertation. 


CHAPTER  II 


LITERATURE  REVIEW 


As  a  specific  example,  this  research  considers  the  PACE  flexible  multibody 
configuration  shown  in  Figure  1.  We  are  interested  in  the  “forward”  problem  of 
mathematical  modeling  such  systems,  and  the  “inverse”  problem  of  estimating  model 
parameters  from  response  measurements. 

Kwak  and  Meirovitch^  employ  Lagrange’s  equations  in  terms  of  general 
quasicoordinates  to  develop  models  for  multibody  systems.  In  this  approach,  the  equations 
of  motion  for  each  individual  flexible  body  are  derived  in  terms  of  quasicoordinates.  A 
recursive  kinematic  description  of  the  velocity  vector  of  a  point  in  each  body  in  a 
multibody  chain  is  expressed  in  terms  of  the  velocity  vector  of  the  preceding  body  in  the 
chain.  The  set  of  equations  are  assembled  globally  using  the  kinematical  relationships  and 
the  redundant  coordinates  are  eliminated.  The  resulting  equations  consist  of  nonlinear 
ordinary  differential  equations  describing  translations  and  rotations  and  partial  differential 
equations  describing  elastic  motions.  The  elastic  displacements  are  discretized  using  a 
finite  element  method  to  approximate  the  solution  of  partial  differential  equations  using  a 
finite  system  of  ordinary  differential  equations.  Kwak  and  Meirovitch  developed  a 
perturbation  approach  to  control  this  high  order  nonlinear  system.  The  perturbation 
approach  assumes  the  rigid  body  motions  are  large  compared  to  the  elastic  motions 
allowing  the  problem  to  be  separated  into  two  parts.  The  problem  is  divided  into  a  low 
dimensional  set  of  nonlinear  zero  order  equations  for  rigid  body  motion  and  first  order 
linear  equations  of  high  dimensionality  for  elastic  motions. 

A  similar  structure  to  PACE  is  studied  by  Meirovitch  and  Lim^.  It  consists  of  two 
hinge  connected  flexible  arms  with  one  end  mounted  on  a  rigid  platform  and  the  other  rigid 
end  holding  a  payload.  The  equations  of  motion  for  this  system  are  derived  using  the 
standard  Lagrange’s  equations.  The  partial  differential  equations  describing  the  elastic 


motions  are  discretized  using  a  Ritz  type  approach  using  shape  functions  called  quasi¬ 
comparison  functions.  Quasi-comparison  functions  are  linear  combinations  of  admissible 
functions  and  in  this  case  are  linear  combinations  of  clamped-free  and  clamped-clamped 
shape  functions.  Essentially  the  same  perturbation  approach  is  used  to  separate  the 
problem  into  a  zero  order  problem  for  rigid  body  maneuvering  and  a  first  order  problem  for 
the  control  of  elastic  motion.  The  elastic  motions  are  assumed  to  be  small  in  this  study  as 
well. 


The  large  deformation  of  beams  has  been  studied  extensively  by  D.H.  Hodges  for 
use  in  rotorcraft  applications.  In  one  study,  Hinnant  and  Hodges'*  modeled  a  constitutive 
nonlinear  effect  on  static  deformations.  They  analyzed  actual  experimental  data  for  a 
cantilever  beam  undergoing  large  deformations.  One  part  of  the  experiment  consisted  of 
measuring  the  static  deformation  of  a  uniform  cantilever  beam  with  a  mass  attached  to  the 
tip.  The  static  behavior  of  the  beam  is  sensitive  to  the  value  of  the  stiffness  coefficient  and 
the  classical  linear  model  gave  incorrect  results.  When  the  beam  deformations  are  large, 
the  linear  theory  is  typically  too  soft.  To  account  for  material  nonlinear  effects,  a  simple 
nonlinear  planar  elastic  model  is  used  which  introduced  a  large  deformation  material 
nonlinearity  coefficient  a.  The  strain  energy  with  this  model  is: 

U  1'^  0) 

where  k  is  the  planar  curvature  of  the  beam  and  h  and  I4  are  area  moments  of  inertia  for 
flatwise  transverse  deflection: 


In  another  study  published  by  Hodges^  the  curvature  of  beams  undergoing  large 
deformations  is  examined.  He  developed  a  formulation  which  accounted  for  the  effective 
shortening  of  a  beam  due  to  transverse  deflections.  Figure  2  shows  the  foreshortening 
effect  where  the  distance  x  is  to  a  point  on  the  beam  before  deformation  and  the  distance  x 
is  the  same  point  after  deformation. 


Figure  2.  Foreshortening  due  to  transverse  displacement 


The  expression  for  curvature  he  derived  is: 


d^w 


(4) 


This  formula  is  exact  for  an  elastic  beam  imdergoing  a  planar  deflection  with  small  axial 
strain  compared  to  unity  and  is  valid  for  any  size  deformation. 


Monasa  and  Lewis^  also  studied  the  large  deformation  behavior  of  beams.  They 
included  the  geometric  nonlinear  effect  due  to  the  curvature  of  the  beam  but  used  a 
different  expression  for  curvature  which  is  frequently  found  in  some  calculus  texts  : 


The  coefficient  of  the  nonlinear  term  in  Eq.  (5)  is  of  opposite  sign  and  the  denominator  is 
of  higher  order  than  die  terms  in  Eq.  (4).  This  difference  is  due  to  the  fact  that  in  the 
formulation  of  Eq.  (5),  jc  is  defined  as  the  distance  along  the  deformed  curve  whereas  in  Eq. 
(4),  X  is  measured  along  the  original  position  of  the  beam.  The  curvature  formula,  Eq.  (5), 
must  be  corrected  when  applied  to  the  physical  beam  curvature  for  the  shortened  value  of  x 
due  to  the  deflection  w.  The  authors  used  an  iterative  process  to  adjust  the  deflection  w 
and  the  u  change  in  x.  First  they  assumed  small  loads  and  no  shortening  and  solved  for  the 
deflection  w  using  Eq.  (5)  and  the  moment  curvature  relationship  and  corrected  the  change 
u  using  the  arc  length  formula.  The  calculations  are  repeated  until  the  solution  converged 
on  u.  Then  the  loads  are  gradually  increased  and  the  iterative  process  repeated  for  each 
step  imtil  the  final  load  is  reached.  The  curvature  formula  caused  a  further  problem  during 
the  Runge-Kutta  integration  of  Eq.  (5)  and  the  moment  curvature  expression.  The  limits  of 
integration  over  the  beam,  x/  and  Xf+j,  are  functions  of  the  beam  deformation  and 
unknown. 

Reddy  and  Singh®  used  the  same  definition  of  curvature  (Eq.  5)  in  their  study  of 
large  deflection  and  free  vibration  behaviors  of  elastic  beams.  It  does  not  appear  they 
accounted  for  the  foreshortened  value  of  x  in  their  solution.  They  compared  two  finite 
element  solutions  of  the  differential  equations  of  motion  and  solved  a  number  of  beam 
problems  with  various  boundary  conditions.  They  used  a  conventional  finite  element 
method  which  used  cubic  interpolating  functions  to  approximate  axial  displacement, 
transverse  displacement  and  slope.  They  also  use  a  finite  element  method  called  a  mixed 
method  which  is  based  on  including  the  bending  moment  with  the  deflections  u  and  w  as 
the  dependent  variables.  By  using  the  bending  moment  as  a  variable,  the  order  of  the 


differential  equation  is  reduced  and  therefore  linear  polynomials  could  be  used  to 
approximate  the  transverse  deflection  instead  of  cubic  polynomials.  They  compared  the 
results  of  their  methods  to  several  analytic  and  numeric  solutions  of  single  beams  with 
various  boundary  conditions. 

Epstein  and  Murray^  defined  the  deformation  of  a  beam  using  two  variables,  axial 
displacement  u  and  transverse  displacement  w.  They  derived  the  following  expression  for 
the  physical  curvature  of  the  beam  which  is  used  for  large  deformations  and  axial  strains: 


f  du^ 

f  > 

d  w 

^dw^ 

^dx^  j 

^dx  j 

^dx^  j 

^  (6) 


They  developed  the  equations  of  motion  for  the  beam  deformation  using  the  principle  of 
virtual  work.  They  discretized  the  displacements  u  and  w  using  a  finite  element  approach 
with  cubic  interpolating  functions.  They  verified  their  formulation  worked  for  large 
deformations  by  computing  the  displacements  from  a  thin  cantilever  strip  bent  by  a 
moment  into  a  complete  circle  and  comparing  the  results  to  the  exact  solution.  They  only 
developed  the  equations  for  the  static  deformation  of  the  beam. 

Several  methods  for  solving  static  and  dynamic  responses  with  large  deformations 
of  single  beams  has  been  discussed  with  various  ways  of  accounting  for  geometric 
(foreshortening)  and  constitutive  nonlinear  effects.  The  modeling  of  two  beam  systems 
have  typically  not  included  these  effects.  A  new  approach  is  developed  in  this  dissertation 
to  model  large  deformation  and  nonlinear  effects  on  a  flexible  two  beam  system  and  this 
model  is  used  to  recover  some  of  the  beam  parameters.  This  methodology  can  be 
generalized  for  use  in  other  studies  based  upon  these  developments. 

First  an  alternative  approach  for  solving  large  static  deformations  of  beams  is 
developed.  This  approach  accounts  for  the  foreshortening  of  beams  in  a  manner  consistent 


with  the  exact  curvature  relationships  of  Hodges  Eq.  (4)  and  Epstein  and  Murray  Eq.  (6). 
However,  this  approach  leads  to  a  new  nonlinear  differential  equation  model  and  an 
associated  method  for  solving  the  resulting  two  point  boundary  value  problem.  Next  a 
nonlinear  finite  element  model  is  developed  for  solving  the  free  vibration  and  dynamic 
response  of  the  beam  system.  The  beam  parameters  affecting  the  deformation  shapes  are 
recovered  separately  for  the  static,  free  response  and  dynamic  response  of  the  system. 


CHAPTER  III 


NONLINEAR  STATIC  ANALYSIS 

The  static  analysis  of  beam  deformation  would  be  the  first  logical  step  in  isolating 
beam  properties.  Beam  parameters  such  as  the  stiffness  coefficient  and  parameterization  of 
the  material  nonlinear  effects  can  be  recovered.  A  mathematical  model  of  the  system  must 
be  developed  before  an  algorithm  can  be  posed  for  recovering  the  beam  parameters.  For 
the  static  analysis,  a  new  approach  for  solving  large  beam  deformations  is  developed. 

3.1  Model  Development  -  Arc  Length  Approach 

A  novel  coordinate  system  is  used  to  develop  the  equation  of  motion  for  an  elastic 
beam  undergoing  transverse  vibration.  This  approach  is  designed  to  accommodate 
geometric  nonlinearities  associated  with  large  deformations.  In  the  conventional  analysis 
the  deformation  is  expressed  by  defining  the  transverse  displacement  at  a  point  x  by  w(x,t) 
and  u(x,t)  shown  in  Figure  2.  In  the  approach  of  this  dissertation,  arc  length  5  is  used  as  the 
independent  distance  coordinate  and  the  slope  angle  is  used  as  the  instantaneous 

generalized  coordinate  shown  in  Figure  3. 


Figure  3.  Transverse  displacement  in  arc  length  coordinates 


The  arc  length  is  the  distance  measured  along  the  instantaneous  deformation  curve 
from  the  origin  and  ^(s,t)  is  the  angle  of  the  local  tangent  vector  e, .  The  vector  r(s,t)  is 

the  inertial  position  of  a  typical  mass  element.  A  tangent  and  normal  vector  is  defined  at 
each  point  along  the  curve.  Figure  4  shows  the  free  body  diagram  for  a  beam  element  of 
length  ds . 


M(s,t)  . 


Q(s.t) 


ds 


Figure  4.  Element  freebody  diagram 


The  normal  force  equation  of  motion  for  this  element  is 


Q{s,t)  +  — ^ — ds 


Q(s,t)  =  m{s)ds 


^^r{s,t)  ^ 


2  -n 


(7) 


and  the  moment  equation  of  motion  is 


^M(s,t) 

Mis,t)  + - — — ds 


-  M{s,t)  + 


Q{s,t)  + - ^ — ds 


l/j  =  I  {s') - — ds 


(8) 


In  the  statics  case,  Eq.  (7)  can  be  reduced  to 
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^Q{s,t) 


=  0 


(9) 


and  by  ignoring  second  order  terms  and  rotatory  inertia  I(s),  Eq.  (8)  can  be  reduced  to  the 
familiar  result 


Qi.s,t)  =  - 


(10) 


The  bending  moment  can  be  related  to  the  curvature  of  the  beam  by  the  Euler  Bernoulli 
constitutive  model 


M(s,t)  =  EI(s) 


where  is  the  local  radius  of  curvature 


(11) 


From  calculus^,  the  curvature  at  a  point  on  a  curve  is  the  derivative  of  the  tangent  angle 
^(s,t)  with  respect  to  the  arc  length  s.  For  the  statics  case,  the  angle  is  denoted  as  T^(5)  and 
the  curvature  is  the  reciprocal  of 


■1  ^  (12) 
^  ds 

In  the  developments  that  follow,  no  small  deformation  approximations  are  used.  This 
exact  curvature  formula  is  combined  with  Eq.  (9)  -  Eq.  (11)  to  yield  the  elegantly  simple, 
•  exact  differential  equation  for  the  static  deformation  of  an  Euler  Bernoulli  beam 


d_ 

ds 


El - ^ 


V. 


ds^ 


=  0 


J 


(13) 


The  (x,y)  shape  of  the  beam  is  obtained  by  solving  the  exact  geometric  differential 
equations 
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^  =  sm(>PW) 

as 

^  =  cos('P(.)) 


(14) 

(15) 


The  choice  of  'P  as  the  dependent  variable  and  s  as  the  independent  variable  has 
two  advantages.  The  first  is  the  "foreshortening"  of  the  beam  is  automatically  and  exactly 
accounted  for.  This  allows  large  deformations  without  having  to  introduce  approximate 
corrections  in  the  length  of  the  beam,  and  furthermore,  simplifies  the  numerical  integration 
since  the  arc  length  is  varied  between  zero  and  the  length  of  the  beam  regardless  of  the 
projected  x  direction  length.  The  second  advantage  is  that  it  allows  the  Euler  Bernoulli 
based  development  to  extend  over  much  larger  geometrically  nonlinear  deformations  since 
the  curvature  of  the  beam  is  not  approximated. 

The  behavior  of  any  beam  will  exhibit  some  material  nonlinearity.  A 
representative  nonlinear  term  is  included  in  the  governing  differential  equation  to  illustrate 
how  such  constitutive  nonlinearities  are  accounted  for.  The  form  of  the  nonlinear  term 
adopted  is  consistent  with  the  model  proposed  by  Hinnant  and  Hodges'*.  The  bending 
moment  that  corresponds  to  the  potential  energy  term  in  Eq.  (1)  is 


A/  =  £7— +  q 

as  \  as  j 


ds 


(16) 


Using  this  definition  of  the  moment  M  and  Eq.  (9)  and  Eq.  (10),  the  static  equilibrium 
equation  becomes 


El 


— ^+6C— 
ds^  ds 


ds^ 


+3d 


vrfj  ds^ 


(17) 


This  cubic  nonlinear  effect  corresponds  to  a  Tiard'  (C>0)  or  'soft'  (C<0)  nonlinear 
stiffening  with  large  deformations. 

3.1.1  Model  Validation 

During  the  1970’s,  Dowell,  Traybar  and  Hodges^°  performed  a  series  of 
experiments  at  Princeton  University  on  the  large  deformation  of  a  cantilevered  beam  with  a 
mass  attached  to  the  tip.  They  were  primarily  interested  in  helicopter  rotor  blade  stability, 
and  carried  out  several  experiments  on  a  cantilever  beam  by  rotating  the  beam  root  at 
various  angles  with  respect  to  the  beam  principal  axis  and  varying  the  tip  weight.  They 
compared  their  results  to  the  Hodges  and  Dowell"  nonlinear  theory  of  rotor  blade 
dynamics,  and  presented  some  results  for  pitch  angles  of  90  degrees  which  corresponds  to 
planar  deflection.  The  arc  length  approach  can  be  checked  against  the  transverse  deflection 
results.  Figure  5  shows  tip  deflections  recovered  using  the  present  arc  length  approach 
along  with  those  predicted  by  the  Hodges  and  Dowell  theory;  and  the  deflections  obtained 
experimentally. 
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Figure  5.  Comparison  of  cantilever  beam  deflections  from  tip  weights 


^  Unfortunately  Dowell,  Traybar  and  Hodges  did  not  report  results  of  beams  at  a  90  degree 

pitch  angle  for  larger  tip  loads,  so  we  cannot  compare  against  their  data  for  larger 
deflections. 

•  However,  published  “exact”  results  are  available  for  large  deformation  of  cantilever 

beams  for  the  special  case  of  a  pure  moment  applied  to  the  free  end.  Figure  6  shows  the 
displacement  profile  of  the  beam  subjected  to  different  moments  from  the  arc  length 
approach  and  from  the  exact  solution’.  These  solutions  allow  the  beam  to  be  deformed  to  a 
A  complete  circle. 


Figure  6.  Deflection  of  cantilever  beam  due  to  tip  moment 


The  arc  length  approach  for  small  and  large  deformations  due  to  end  moments  is  in  very 
good  agreement  both  with  available  independent  experiments  and  with  the  pure  moment 
exact  solutions  from  reference  9. 

3.1.2  Application  to  Two  Beam  Model 

The  arc  length  approach  is  applied  to  a  two  beam  system  similar  to  the  PACE 
structure  in  Figure  1.  The  geometry  of  this  two  beam  system  is  illustrated  in  Figure  7.  The 
PACE  structure  is  simplified  in  that  the  shoulder  hub,  elbow  connection  and  tip  masses  are 
treated  as  point  masses.  In  this  discussion,  the  shoulder  hub  is  locked  and  the  elbow  angle 
is  constrained  to  be  zero  (in  the  'elbow  locked'  configuration).  The  arc  length  of  the  first 


beam  is  measured  from  the  shoulder  hub  and  the  first  beam's  tangent  angle  ^1(5)  is 
measured  relative  to  a  horizontal.  For  the  second  beam,  the  arc  length  is  measured  from 
the  elbow  connection  of  the  two  beams  and  the  tangent  angle  'F2(s)  is  measured  relative  to 
a  vector  tangent  to  the  tip  of  beam  1  (vector  has  angle  relative  to  the  horizontal). 

The  choice  of  these  coordinates  simplifies  the  equation  development  of  beam  2  and  also 
permits  automatic  matching  of  the  displacements  and  slopes  at  the  connections  of  the  two 
beams  (for  a  “locked’  elbow  with  zero  elbow  angle).  This  choice  of  coordinates  allows  the 
Euler  Bernoulli  beam  theory  to  apply  over  large  deformations.  Notice  that  the  angle  ^^2(5) 
is  measured  relative  to  the  tip  of  the  first  beams’  tangent  vector,  whose  elevation  angle  is 
^^^(Ll)  and  is  not  measured  relative  to  the  horizontal. 


Figure  7.  Two  beam  configuration 


The  static  equilibrium  equations  are  derived  for  each  beam  and  differential  equations  for 
the  coordinates  x  andy  are  included  since  they  are  needed  in  the  enforcement  of  a  boundary 
condition: 


Beam  1 


Ell — 7-^ —  +  6C1  - 


ds:^ 


ds 


K  *1^  } 


+  3C, 


d^lisQY  d^^iisQ  _ 

dsi  J  dsj^ 


0  (18a) 


4^1  (■^1) 
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=  5111(^1(51)) 


(18b) 


^l(^l) 

dsi 


=  COS(^i(5i)) 


■(18c) 


Beam  2 


d^W2(s2)  ,  _  ^^2(^2) 
El  7 - : — 5 —  +  0C2 


ds2- 


ds2 


^d^^2M 
ds2^ 


^2 


+  3Q 


r^2(^2)r^  ^2(-^2)_q 

V  ^2  j  ds2^ 


(19a) 


dyiM 

ds2 


=  sin(^2(-52)) 


(19b) 


=  C0S('F2  (^2  ))  (  ^ 

dS2 

There  are  a  total  of  ten  boundary  conditions  that  the  solution  of  the  above  tenth 
order  system  of  differential  equations  must  satisfy.  Six  boundary  conditions  apply  at  the 
origin  of  each  beam  (5  =  0)  where  the  angle  ^(5)  and  the  coordinates  x(0),  XO)  must  equal 
zero.  Two  boundary  conditions  apply  at  the  elbow  connection  of  the  beams  where  the 
bending  moments  must  be  equal  Eq.  (20)  and  the  shear  forces  must  be  equal  Eq.  (21). 
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(21) 


Two  final  boundary  conditions  must  be  enforced  at  the  tip  of  the  beam  system  depending 
on  how  the  boundary  forces  are  applied.  The  wire  constraint  experiment  provides  the  final 
two  boundary  conditions  needed  to  solve  for  the  beam  deformation  shape. 

3.1.3  Wire  Constraint  Experiment 


A  nonlinear  static  deformation  experiment  has  been  conceived  to  simplify  testing  of 
the  two  beam  system  and  allow  recovery  of  beam  properties  which  contribute  to  the 
potential  energy  of  the  system.  The  idea  of  the  experiment  is  to  attach  a  wire  to  the  tip  of 
the  beam  system  and  measure  the  force  required  to  deform  the  beam.  A  series  of 
experiments  with  different  measured  tensions  in  the  wire  would  yield  a  family  of  measured 
deformation  shapes  from  which  beam  parameters  such  as  the  stiffness  coefficients  could  be 
recovered.  In  most  cases,  only  the  tip  and  elbow  deflections  need  to  be  measured  for 
estimating  the  static  parameters.  The  use  of  a  wire  with  negligible  bending  stiffness  has 
several  advantages.  From  the  beam  standpoint,  the  wire  provides  a  point  force  but  no 
moment  which  acts  on  the  tip  of  the  beam.  Along  the  wire,  this  one  dimensional  force  can 
be  measured  anywhere  using  a  load  cell.  The  unit  that  holds  the  load  cell  and  the  end  of 


the  wire  would  be  attached  to  the  tip  of  the  beam  and  to  the  horizontal  table  using  suction 
cups  or  electromagnets.  This  would  also  allow  the  beam  to  be  released  quickly  without 
interfering  with  the  free  response  of  the  beam.  Once  release  occurs,  the  wire  could  be 
removed  quickly  allowing  a  free  vibration  of  the  two  beam  system  and  the  recovery  of 
some  dynamic  properties.  Thus,  the  release  mechanism  serves  two  functions:  (i) 
measurement  of  static  deformation  forces  and  (ii)  a  release  mechanism  for  free  dynamic 
response  experiments.  A  prototype  mechanism  which  measures  the  static  tip  force  and 
permits  easy  release  of  the  wire  has  been  designed  and  tested  by  I.  Romero  . 

Radial  Wire  Constraint 

One  proposed  experiment  uses  the  wire  to  connect  the  tip  of  the  beam  to  the 
shoulder  hub.  The  wire  would  attach  to  a  slip  collar  arotmd  the  shoulder  hub  and  only  the 
magnitude  of  the  force  would  be  measured.  A  series  of  experiments  with  different  wire 
lengths  and  corresponding  required  forces  would  be  used  to  recover  the  parameters.  For 
this  static  experiment,  the  two  beam  system  must  satisfy  the  geometric  constraint  due  to  the 
fixed  length  of  the  wire.  This  fixed  wire  length  has  an  associated  radial  constraint  force 
(wire  tension)  required  for  the  beam  deformation.  Figure  8  illustrates  the  geometry  of  the 
wire  length  constraint  for  the  special  case  that  the  opposite  end  of  the  wire  with  the 
measurement  load  cell  is  attached  to  the  shoulder  joint.  This  “bow”  setup  is  but  one  of  an 
infinite  number  of  choices . 
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The  following  equations  define  the  geometric  constraint 


where 

Xt  =  +  Xf2  cosC'FiCI-i))  ->'/2  sinCTiCIi))  (22b) 

and 

yt  =  yt\  +  sin('Fi(Z,i))  +  yti  cos('Pi(Ii))  (22c) 

The  radial  wire  provides  the  final  two  boundary  conditions;  the  geometric  constraint 
Eq.  (22a)  -  Eq.  (22c)  must  be  satisfied  and  the  bending  moment  at  the  tip  of  the  second 


beam  must  be  zero. 


generated  by  simply  deforming  the  beam  with  different  wire  tensions.  Applying  a  directed 
force  would  be  more  difficult  to  implement  experimentally  since  the  shoulder  hub  must  be 
constrained  and  the  alignment  of  the  force  would  have  to  be  done  carefully.  For  this 
experiment,  the  final  two  boundary  conditions  would  require  the  shear  force  in  the  wire  to 
match  the  shear  force  in  the  beam 


=F^sin('P2(i^2)  +  ^l(A))  =  -E^2— ^r^  +  3C2 


^2(^2) 
ds') 


j2 


ds2^ 


(27) 


and  the  bending  moment  at  the  tip  is  zero,  Eq.  (23). 

A 

3.1.4  Boundary  V alue  Problem 

Satisfying  the  set  of  equilibrium  equations  Eq.  (18a)  -  Eq.  (19c)  is  a  nonlinear  two- 
point  boundary  value  problem  where  some  initial  conditions  are  unknown.  An  algorithm 
for  the  solution  of  this  boundary  value  problem  is  established  as  follows.  The  first  beam 
differential  equations  Eq.  (18a)  -  Eq.  (18c)  are  integrated  using  the  Runge  Kutta  technique 
along  the  length  of  the  beam,  then  the  second  beam  equations  Eq.  (19a)  -  Eq.  (19c)  are 
integrated.  A  second  order  Newton's  method  is  used  to  iteratively  adjust  approximate 
initial  conditions  to  match  the  boundary  conditions.  The  initial  conditions  are  summarized 
below 

Initial  Conditions  -  Beam  1 

'Pl(0)  =  0  (28a) 


^l(O) 

dsi 


=  unknown 


(28b) 
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• 

rf'ri2(0) 

.  -  unknown 

dsi^ 

(28c) 

Tl(0)  =  0 

(28d) 

• 

xi(0)  =  0 

(28e) 

• 

Initial  Conditions  -  Beam  2 

'P2{0)  =  0 

(29a) 

• 

dVjiO)  . 

-  unknown 
dS2 

(29b) 

• 

=  unknown 

ds2^ 

(29c) 

y2(0)  =  0 

(29d) 

• 

X2(0)  =  0 

(29e) 

In  order  to  determine  the  unknown  initial  conditions  (^i'(O)  ,'Fi"(0),'F2'(0)  using 

Newton's  method,  the  derivatives  of  the  boundary  conditions  [Eq.  (20),  Eq.  (21),  Eq  (22), 
Eq.  (23)]  with  respect  to  these  parameters  are  taken.  The  four  boundary  equations  are 
rewritten  so  that  all  the  terms  are  collected  on  one  side  and  then  the  equations  are 
represented  by  the  variables  Ti,  r2,  Ts,  r4.  The  derivatives  of  these  boundary  conditions 
are  taken  with  respect  to  the  initial  conditions  and  the  results  are  rearranged  and  written  in 
matrix  form  to  solve  for  approximate  imtial  boundary  condition  corrections  as 
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(30) 


The  new  values  of  the  initial  conditions  are  determined  and  the  numerical  solution  of  the 
equilibrium  equations  is  repeated  until  the  unknown  initial  conditions  have  converged.  In 
the  derivatives  of  these  boundary  conditions,  the  sensitivity  of  the  solution  variables,  ('Pj^ 
'P2)  with  respect  to  the  unknown  initial  conditions  are  contained.  For  example,  one  of  the 
terms  needed  in  Eq  (30)  is 
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The  sensitivities  of  the  solution  variables  with  respect  to  the  initial  conditions  (like 
^'2(12)/ ^  2(0))  must  be  numerically  integrated.  This  is  done  using  a  state 
transition  matrix  approach.  One  can  write  Eq.  (18a)  -  Eq.  (18c)  in  first  order  form  using 
the  vector  z(s)  [where  z^s)^  ={'Pi(5),  'PiX^),  'Pr(-s)»  yii^),  ^1(5)}]  and  by  defining  the 

derivative  equations  as  the  vector y(^(s))  so  that 

^  =  ms))  (32) 


By  manipulating  this  equation,  the  following  differential  equation  is  developed  which 
governs  the  spatial  evolution  of  the  "transition  matrix"  (partial  derivatives  of  z(s)  with 
respect  to  z(0)): 


d 

<^(5) 

ds 

_  ^2(0)  _ 

_  ^2(0) _ 

(33) 


The  initial  (s=0)  condition  for  this  matrix  of  equations  is  the  identity  matrix  and  these 
equations  can  be  numerically  integrated  simultaneous  with  Eq.  (18a)  -  Eq.  (18c)  over  the 
arc  length  to  yield  the  sensitivities  to  'Fi'(O)  and  T'i"(0).  The  same  approach  can  be  taken 
for  the  second  beam  equations  [Eq.  (19a)  -  Eq.  (19c)]  using  the  analogous  state  vector  w(s) 
and  defining  the  derivative  equations  as  the  vector  g(w(s)) 

Manipulation  of  this  equation  yields 


d 

^(s) 

ds 

_^w(0)_ 

Sv{s) 

_  ^(0)_ 

(35) 


The  initial  condition  of  this  set  of  equations  is  the  identity  matrix  and  the  equations  are 
numerically  integrated  along  the  beam  length.  Values  for  the  terms  ^2(^2)/ ^z(O)  and 
Sv(0)  are  used  in  the  derivatives  of  the  boundary  equations  (in  terms  like 
Eq.  (31)),  and  the  changes  to  the  unknown  initial  conditions  are  calculated. 


3.1.5  Simulations 


Radial  Wire  Constraint 

The  deformation  of  the  two  beam  system  with  beam  properties  from  the  PACE 
model  is  studied  for  various  lengths  of  the  wire  constraint  and  values  of  the  nonlinear 
coefficients  Cj  and  C2.  The  iterative  solution  process  outlined  in  the  previous  section 
converges  very  rapidly  and  without  any  numerical  singularities  for  wire  lengths  ranging 
from  98%  to  42%  of  the  total  beams  length.  Obviously  the  traditional  Euler  Bernoulli 
beam  theory  would  not  apply  to  some  of  these  large  deformations  but  it  does  show  the 
robustness  of  the  solution  technique.  Figure  9  shows  the  deformation  of  four  cases  with 
different  wire  constraints. 


Figure  9.  Effect  of  wire  constraint  length 


The  force  in  the  wire  increases  dramatically  with  successively  smaller  wire  lengths.  Figure 
10  shows  the  change  in  deformation  shape  for  a  given  length  of  wire  with  increases  in  the 
nonlinear  coefficients.  These  results  are  for  positive  C's  which  correspond  to  beam 
stiffening. 


The  first  increase  in  the  coefficients  of  the  nonlinear  terms  results  in  the  largest  change  in 
the  deformation  shape  of  the  system.  Larger  increases  in  the  coefficients  slightly  change 
the  deformation  shape  but  the  effects  of  the  increases  are  mainly  felt  in  the  wire  force 


calculations  which  involves  higher  derivatives  of  the  angle  Table  2  shows  the  trend 

of  larger  increases  in  the  wire  force  and  smaller  effects  on  the  systems  tip  location. 


Table  2.  Material  nonlinearity  results 

Coefficient  Cj,  C2 

Force  in  wire 

Tip  X  location 

Tip  y  location 

0  (No  nonlinearity) 

33.3 

.506 

1.197 

1,1 

37.6 

.573 

1.167 

5,5 

53.3 

.643 

1.130 

10, 10 

72.1 

.668 

1.115 

20, 20 

109.3 

.686 

1.104 

30,  30 

145.3 

.693 

1.100 

Directed  Force 

In  this  experiment,  the  force  on  the  tip  of  the  system  is  directed  perpendicular  to  the 
undeformed  axis  (the  y  direction  shown  in  Figure  8).  This  force  will  be  referred  to  as  a 
transverse  force.  The  deformation  of  the  two  beam  system  with  beam  properties  from  the 
PACE  model  is  studied  for  different  values  of  wire  forces.  To  be  consistent  with  the  initial 
deformation  shape  used  in  the  dynamic  analysis,  the  nonlinear  coefficients  are  set  to  zero. 
This  removes  the  constitutive  nonlinearity  from  the  model  but  leaves  the  geometric 
nonlinearity  due  to  the  trigonometric  functions  involved  in  the  solution.  The  deformation 
shapes  for  four  wire  tensions  are  shown  in  Figure  11. 
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Figure  11.  Deformation  shapes  for  series  of  transverse  forces 


The  solution  of  the  boimdary  value  problem  quickly  recovers  the  unknown  initial 
conditions  and  the  deformation  shapes  for  the  different  force  boundary  conditions. 

3.2  Static  Parameter  Estimation 

If  the  beam  properties  are  uncertain,  the  parameters  affecting  the  static  deformation 
shape  can  be  recovered  using  measurements  from  two  or  more  experiments.  The  effective 
stiffness  and  nonlinear  coefficients  for  each  beam  are  some  of  the  parameters  that  may  not 
be  well  known. 


3.2.1  Least  Squares  Estimation 


A  least  squares  analysis  is  used  to  recover  the  best  fit  parameters  for  the  arc  length 
models  represented  by  Eqs.  (18a)  -  (18c)  and  Eqs.  (19a)  -  (19c).  This  is  a  nonlinear 
estimation  problem  and  the  least  square  differential  correction  technique^^  is  used  to 
successively  approximate  the  parameters.  The  error  between  the  measured  values  and 
those  predicted  by  the  model  using  the  updated  parameters  are  minimized 

^  ~  '^meas  ~  ^ model  ^Pn  +l)  0^) 

The  parameter  vector  p  represents  the  variables  to  be  recovered.  Measurements  at  different 
time  intervals  are  represented  by  the  vector  .  The  model  values  cannot  be  written  as 
explicit  functions  of  the  parameters  since  determining  the  model  values  involves  the 
numerical  solution  of  the  differential  arc  length  equations.  The  values  from  the  nonlinear 
model  for  a  given  set  of  parameters  can  be  calculated  and  a  new  set  of  parameters  can  be 
approximated  using  a  first  order  Taylor  series  expansion 


YmodeliPn+l)  =  ^model(Pn)  +  ^(pn+l-Pn)  .  ^here  Q  = 


(37) 


This  is  substituted  into  the  error  expression  giving 
^  “  ''^meas  ~  ^ model  ^Pn )  ~  +1  Pn^ 


(38) 


The  sum  square  of  the  error  is  minimized  and  the  new  approximation  for  the  parameter;?  is 
solved 


Pn+l  ~  Pn  ^  ^ meets  model  ^Pn^) 


(39) 
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Each  column  of  the  sensitivity  matrix  Q  represents  the  derivative  of  the  corresponding 
measurements  with  respect  to  one  parameter.  The  rows  of  the  matrix  represent  the  values 
of  the  gradients  at  each  measurement  time. 

This  estimation  procedure  is  complicated  by  the  fact  that  the  governing  equations 
are  nonlinear  differential  equations  and  further  complicated  by  the  fact  any  trial  solution 
requires  solution  of  a  two  point  bomdary  value  problem.  In  order  to  establish  the  needed 
sensitivities,  we  introduce  a  novel  method  for  differentiating  an  implicit  function  . 

Implicit  Function  Theorem 

.  The  differential  equations  must  first  be  solved  for  the  initial  values  of  T  for  a.  given 
set  of  parameters  p  and  then  repeated  for  updated  values  of  the  parameters  until  the  method 
converges.  A  more  subtle  question  arises,  what  is  the  sensitivity  of  the  solution  ('Pi'(O), 
T^l"(0),  'F2'(0),  'P2"(0)  )  of  the  boundary  value  problem  with  respect  to  a  parameter  pi 
Therefore  the  boundary  value  problem  must  be  solved  for  each  iteration  of  the  parameters. 
In  order  to  calculate  the  changes  in  the  parameters  using  the  least  squares  method,  the 
derivatives  of  the  measurements  (such  as  xf^  with  respect  to  the  parameters  must  be 
calculated.  To  fill  this  derivative  matrix,  the  equations  defining  the  measurements 
(Eq.  (22b)  for  Xf )  must  be  differentiated  with  respect  to  the  parameters  p.  These  derivative 
equations  contain  terms  such  as  dx\{L])l  dp,  d^\{L\)l  dp  which  must  be  determined  and 
integrated  along  the  arc  length.  Again  a  state  transition  approach  is  taken  but  in  this  case 
we  are  concerned  with  the  model  parameters  p  instead  of  the  initial  conditions  z(0).  The 
governing  differential  equations  for  the  first  beam  Eq.  (18a)  -  Eq.  (18c)  can  be  written 

^(z(s))=/(z,p)  (40) 

and  by  manipulation,  we  can  derive  a  differential  equation  for  the  needed  partial 
derivatives: 
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•  A. 

>z(5)' 

r  ^/i 

ds 

^  . 

(41) 


# 
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Similarly  the  governing  equations  for  the  second  beam  Eq.  (19a)  -  Eq.  (19c)  can  be  written 


ds 


{w{s))  =  g{w,p) 


and 


d 

^w(s) 

ds 

^P 

yp\ 

4? 


(42) 


(43) 


To  integrate  these  matrix  differential  equations  to  compute  the  derivatives  {dx\{L\)ldp, 
d^\(L\)ldp  ,  etc),  the  initial  conditions  dz{(i)ldp  and  dw{0)ldp  must  first  be  determined. 
The  main  difficulty  lies  in  determining  d^\{())ldp,  dV\"{(!i)ldp,  d'V'^{Q)ldp  and  d'^-^XOydp 
which  are  the  initial  conditions  of  these  sensitivity  matrices  in  Eq.  (41)  {dz{(S)ldp)  and  Eq. 
(43)  {dw{0)ldp).  The  initial  values  'Pi'(O),  'Pi"(0),  'F2'(0)  and  'F2"(0)  must  first  be 
determined  to  match  the  boundary  conditions  as  specified  in  section  3.1.4  (using  the 
differential  equations  Eq.  (18a)  -  Eq(19c)).  Then  the  initial  conditions  for  the  state 
transition  matrices  Eq.  (41)  and  Eq.(43)  can  be  determined.  To  find  the  initial  conditions 
for  d^xiOydp,  a'Fi"(0)/S/?,  d^i(Qi)ldp  and  d^2\^)ldp  wq  will  apply  the  implicit  function 
theorem  to  the  boundary  condition  equations.  The  boimdary  condition  equations  [Eq.  (20), 
Eq.  (21),  Eq.  (22)  and  Eq.  (23)]  can  be  written  as  functions  of  the  initial  conditions  and  the 
parameters/? 

ri('Pi’(0),'Pi"(0),'P2'(0),^2’'(0),£fi,ii:i,^/2.^2)  =  o  (44) 

r2('Pi'(0),Yl”(0),'F2'(0),'P2''(0).^A.^l>^^2.^2)  =  0  (45) 


r3(Yl'(0),Ti"(0),'P2'(0),'P2”(0),^/l,^l,i?^2.^2)  =  0 


(46) 
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r4('Pi-(0),'Pi"(0),'P2’(0X^2"(0X^A.^b^^2.^2)  =  0  (47) 

We  can  treat  the  initial  conditions  Ti'(O),  'Fi"(0),  'P2'(0).  'P2"(0)  as  the  dependent 
®  variables  and  the  parameters  p  as  the  independent  variables.  The  implicit  function  dieorem 

is  used  to  determine  derivatives  of  these  functions.  This  is  an  unusual  application  of  the 
theorem  since  the  variables  in  the  boundary  condition  equations  are  obtained  from  a 
solution  of  a  two  point  boundary  value  differential  equation.  Using  the  derivative  of  an 
^  implicit  function^,  the  partial  derivatives  of  the  functions  with  respect  to  each  parameter  pi 

is 


. 

^.'(0)  . 

^i"(0) 

^2'(0) 

^2"(0)  . 

• 

4^,-  "  ^,'(0) 

47,  ^2'(0) 
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(48) 
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(51) 


The  implicit  and  explicit  derivatives  are  computed  and  inserted  into  this  set  of  equations 
which  are  used  to  solve  for  <^{(0)/  . 

These  are  needed  along  with  other  initial  conditions  to  integrate  the  sensitivity  matrices 
Eq.  (41)  and  Eq.  (43)  over  the  arc  lengths  5.  After  this  integration,  the  recovered 
derivatives  such  as  dxi{L{)ldEIi  are  used  in  computing  the  measurement  equation 
derivatives. 

The  computation  of  the  derivatives  in  the  sensitivity  matrix  can  now  be  carried  out 
knowing  the  partials  of  the  first  and  second  beam  variables  from  the  process  outlined 
above.  With  this  derivative  matrix  we  can  compute  the  needed  corrections  to  the 
parameters  p  using  Eq.  (39).  This  whole  process  must  be  repeated  until  the  parameters  p 
are  recovered  which  yield  the  measurement  values  for  each  test.  As  will  be  evident  in  the 
numerical  studies,  the  domain  of  practical  convergence  is  usually  large,  tolerating  up  to 
40%  errors  in  the  starting  estimates. 


Radial  Wire  Constraint 

The  radial  wire  experiment  permits  us  to  estimate  the  effective  stiffness  and 
nonlinear  coefficients  {p  =  Eli,  ■^^2’  ^2)  ^sing  as  few  as  two  experiments  with  different 
wire  lengths.  There  are  two  independent  measurements  taken  during  each  test.  In  this 
development,  the  force  in  the  constraining  wire  and  the  global  x  location  of  the  tip  of  the 
second  beam  are  assumed  to  be  measured.  Any  x  or  y  measurement  would  also  work,  but 
Xf  is  chosen  for  experimental  simplicity.  For  this  determined  system,  the  sensitivity  matrix 
for  the  two  experiment  case  is 
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M2 

My,^ 
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Ml 

MI 2 

M2 

(52) 


( where:  superscript  ^  denotes  data  from  the  1st  test,  and  ^  denotes  data  from  the  2nd  test) 


The  derivatives  of  the  equations  for  x,  [Eq.  (22b)]  and  theFw  [Eq.(26)]  equations  are  taken 
with  respect  to  the  parameters.  The  partials  of  the  first  and  second  beam  variables 
contained  in  those  equations  are  determined  from  the  process  outlined  in  the  previous 
section.  These  derivatives  are  used  to  populate  the  sensitivity  matrix  and  the  changes  to  the 
beam  parameters  p  are  determined  using  the  least  squares  method. 

Measurement  Errors 

Often  the  measurements  in  Y^eas  made  with  unequal  precision  and 
approximate  weighting  should  be  included  in  the  parameter  updating.  The  reciprocal  of  the 
error  variance  is  the  conventional  “optimal”  choice  for  the  weight  and  is  qualitatively 
correct  since  a  measurement  with  a  small  error  would  have  a  large  weight  and  one  with  a 
large  error  would  have  a  very  small  weight.  The  error  variance,  is  approximated  in 
practice  as  fiie  square  of  the  standard  deviation,  representing  the  theoretical  squared  sums 
of  the  differences  between  a  large  number  of  measurements  and  the  mean  value.  The  least 
squares  function  of  the  residual  error  including  a  positive  definite  weight  matrix  W 
becomes 


*/  =  ~  j^Y^g^  -  ^ model  (Pn  +1 )]  ^ model  iPn  +1 )_ 


(53) 
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where  PF  = 


0 


0 

1 


0 


0 


4 


meas  J 


Substituting  in  for  Ymodel  (Pn  +l)  minimizing  the  function  J,  the  parameter  update 
becomes 


Pn+l  =Pn+  (Cl^WQr^Q^W(Y^,as  -  ^model  iPn)) 


(54) 


The  measurement  errors  are  mapped  through  the  estimation  algorithm  into 
associated  errors  in  the  estimated  parameters.  This  relationship  is  developed  in  reference 
13.  The  covariance  matrix  of  the  parameters  is  related  to  the  variance  matrix  of  the 
measurements  by  the  following 


E,=(Q^fm) 


-’I 


<^pl  ^p\ 


(^pl 


^p#  ^pl  %2  ^p#  ^p2 


%  ^p\  ‘^p# 
^p2  *^p# 


(55) 


The  diagonal  elements  of  the  error  covariance  matrix  gives  an  indication  of  how  accurately 
the  parameters  are  estimated.  Each  of  the  ‘true’  parameters  should  be  typically  contained 
within  three  times  the  standard  deviation  of  the  recovered  parameter.  The  ’s  are 

measures  of  correlation  of  the  estimation  errors. 
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Directed  Force 

The  directed  force  experiment  recovers  the  effective  stiffness  coefficients 
ip  =  Ell,  ^sing  a  number  of  measurements.  In  this  development,  the  global  y  location 

of  the  tip  of  the  second  beam  is  assumed  to  be  measured.  The  y  measurement  is  chosen 
since  the  stiffiiess  coefficients  clearly  have  a  large  affect  on  this  variable.  A  number  of 
measurements  are  used  to  over  determine  the  parameters.  The  sensitivity  matrix  for  this 
problem  becomes 

^yt  ^  ^yt  ^ 

^EIi  ^El2 

^yt  ^  ^yt  ^ 

^EIi  ^El2 

^3  (56) 

^EIi  ^El2 

^yt  ^  ^yt  ^ 

mi  mi 

( where  the  superscript  number  refers  to  the  measurement  number  from  1  to  iV ) 

The  derivatives  of  the  equation  for  y,  [Eq.  (22c)]  are  taken  with  respect  to  both  stiffness 
parameters.  The  partials  of  the  first  and  second  beam  variables  contained  in  those  equations 
are  determined  from  the  implicit  differentiation  process  outlined  previously.  These 
derivatives  are  used  in  the  sensitivity  matrix  and  changes  to  the  beam  parameters  p  are 
determined  using  the  least  squares  method. 

The  parameter  estimation  method  for  the  directed  force  experiment  is  enhanced  to 
include  measurement  errors  and  their  associated  variance  in  the  least  squares  determination 
of  the  parameters.  A  weight  matrix  is  included  in  the  parameter  updating  and  reflects  the 
uncertainties  in  the  measurements.  The  parameters  recovered  have  associated  error 
variances  which  indicates  how  well  the  experiment  determines  the  parameters. 


3.2.2  Results  -  Static  Deformation  Parameters 


The  approach  for  estimation  of  parameters  affecting  the  static  deformation  of  the 
beam  system  is  tested  using  a  computer  simulation  of  measured  data  for  two  types  of 
experiments. 

Radial  Wire  Constraint 

We  first  consider  the  determined  case,  a  radial  wire  experiment  to  determine  the 
static  potential  energy  parameters,  assuming  the  system  measurements  are  accurate.  This 
proposed  experiment  uses  a  wire  to  connect  the  tip  of  the  beam  to  the  shoulder  hub  which 
provides  a  geometric  constraint  on  the  solution.  To  generate  simulated  measurements,  the 
two  beam  forward  solution  using  the  arc  length  approach  is  computed  for  two  different 
wire  lengths  and  a  given  set  of  'true'  parameters.  The  deformation  shapes  of  the  two  wire 
lengths,  Lw=1.4  and  Lw=L3,  are  shown  in  Figure  9.  The  resulting  values  of  Xf  and  are 
used  as  measurements  in  the  inverse  solution  to  see  if  the  program  can  recover  the  true 
values  given  poor  starting  estimates.  An  approximation  of  the  true  model  is  constructed  by 
increasing  the  parameter  errors  from  8%  to  11%  and  is  used  to  start  the  parameter 
estimation  process.  Table  3  summarizes  the  input  and  recovered  values  of  the  parameters. 
The  program  converges  rapidly  to  the  values  in  the  last  column  which  are  virtually 
identical  to  the  'true'  values.  The  calculated  values  of  X(  and  are  matched  exactly  with 
the  input  measured  values. 


1  Tables.  Parameter  recovery  for  wire  constraint 

Parameters 

True  Values 

Input  Values 

Converged  Values 

Ell 

11.413 

10.400 

11.415 

Cl 

5.000 

4.550 

5.000 

El2 

11.275 

10.175 

11.276 

Cl 

20.000 

17.800 

19.991 

The  parameter  recovery  method  based  on  this  series  of  static  measurements  converges 
rapidly  to  the  known  set  of  parameter  values.  The  parameter  C2  is  less  well  observable, 
however,  and  the  true  measurement  case  suggest  larger  deformations  and  redundant 
measurements  will  be  needed  to  accurately  estimate  C2. 

Directed  Force 

This  experiment  recovers  the  effective  stifftiess  coefficients  for  each  beam  when  the 
system  is  subjected  to  a  transverse  force  boundary  condition.  Here  we  ignore  the  terms 
containing  Ci,  C2  and  estimate  only  Eli,  EI2.  The  static  parameters  are  recovered  when  the 
system  is  statically  deformed  in  a  manner  consistent  with  the  initial  conditions  used  jn  the 
dynamic  analysis.  The  measurements  of  y,  are  simulated  by  computing  the  forward 
solution  using  the  arc  length  approach  for  four  values  of  the  transverse  force  F  and  a  given 
set  of  ‘true’  parameters.  This  family  of  static  deformation  shapes  is  shown  in  Figure  11. 
Redundant  measurements  were  used  to  allow  the  overdetermination  of  the  parameters. 
Estimates  of  the  stiffness  coefficients  are  used  to  start  the  parameter  recovery  program  and 
the  results  are  shown  in  Table  4. 


Ta 

ble  4.  Parameter  recovery  for  directed  force 

Parameters 

True  Values 

Input  Values 

Converged  Values 

Ell 

11.413 

10.750 

11.413 

EI2 

11.275 

10.500 

11.275 

The  values  of  the  parameters  are  converged  after  five  iterations  and  are  identical  to  the 
‘true’  values. 

Physical  measurements  of  yt  will  have  some  uncertainties  associated  with  the 
values  and  this  can  be  reflected  in  the  error  variance  of  the  measurements.  To  reflect  this 


in  the  simulated  measurements,  errors  are  introduced  using  a  Gaussian  random  number 
generator  and  an  associated  standard  deviation  a.  For  each  measurement,  an  error  is 
created  from  a  random  number  belonging  to  a  distribution  having  zero  mean  and  the 
specified  standard  deviation  cr.  These  simulated  measurement  errors  are  added  to  the 
‘true’  values  of  the  measurable  quantities  which  are  then  used  in  the  parameter  recovery 
program.  The  weighted  least  squares  function  is  used  to  minimize  the  errors  and  the 
parameters  are  determined  using  Eq.  (54).  The  ‘true’  measurements  of  yt  have  values 
ranging  from  0.276332  to  0.597325.  The  standard  deviations  for  the  errors  in  these 
simulated  measurements  are  input  to  the  program  and  range  from  .0008  to  .0011.  The 
parameter  recovery  program  is  started  with  initial  estimates  of  the  stiffness  coefficients  and 
Table  5  shows  the  recovered  parameters  when  measurement  uncertainties  are  included. 


Table  5.  Parameter  recovery  for  directed  force 
with  measurement  errors 

Ell 

EI2 

True  values 

11.413 

11.275 

Initial  estimate 

10.750 

10.500 

Recovered  values 

11.410 

10.987 

The  stiffness  coefficient  recovered  for  the  first  beam  is  virtually  identical  to  the 
uncorrupted  ‘true’  value  and  has  a  small  variance.  However,  the  stiffness  coefficient 
recovered  for  the  second  beam  is  not  as  well  determined  in  the  presence  of  measurement 
errors  and  differs  by  3%  to  the  ‘true’  value.  The  recovered  parameters  have  associated 
error  covariances  as  defined  by  Eq.  (55).  Small  values  indicate  the  parameters  are  well 
determined.  The  covariance  matrix  for  these  parameters  are 


.004497 

-.047864 


-.047864 

349850 


As  expected,  the  stiffness  coefficient  for  the  second  beam  has  a  large  standard  deviation 
{(j  =  =  .7)  which  indicates  it  is  not  as  well  determined  as  the  first  beams  value.  We 

note  the  converged  estimate  differs  from  the  true  value  of  EI2  by  =  0.28  which  is  well 
within  one  sigma,  so  there  is  good  consistency  between  the  actual  estimation  error  and  the 
uncertainty  estimate.  The  difference  in  the  ability  to  estimate  the  first  beam’s  stiffness 
coefficient  better  than  the  second  beam’s  stiffness  coefficient  is  due  to  the  storage  of 
potential  energy  in  a  cantilever  beam  system.  When  the  system  is  subjected  to  a  tip  force, 
the  portion  of  the  system  closest  to  the  clamped  boundary  condition  deforms  the  most  (has 
the  largest  curvature  and  bending  moment)  and  the  parameters  most  affecting  that  response 
will  most  easily  be  recovered.  The  information  on  the  estimated  parameters  and  their 
associated  covariance  matrix  can  be  “carried  forward”  and  used  in  the  recovery  of 
parameters  affecting  the  dynamic  motion. 


CHAPTER  IV 


NONLINEAR  FREE  VIBRATION  ANALYSIS 

4.1  Model  Development 

A  mathematical  model  for  predicting  the  dynamic  motion  of  the  two  beam  system  is 
developed  for  use  in  estimating  the  parameters  affecting  the  free  vibration  response  of  the 
system.  A  nonlinear  finite  element  model  is  developed  due  to  difficulties  encoxmtered  in 
the  arc  length  approach. 

4.1.1  Arc  Length  Approach 

The  arc  length  approach  can  not  be  easily  extended  to  the  dynamic  analysis  due  to 
the  fact  that  certain  integrals  arise  which  complicate  the  structure  of  the  kinetic  energy  and 
the  equations  which  result.  This  can  be  seen  by  looking  at  the  dynamic  equation  of  motion 


m{s)ds 


+  EI 


=  0 


(57) 


where  the  global  acceleration  is 


r(s,t)  =  x(s,  t)l + Oi 


and  the  unit  normal  is 


= -sin'P(5,r)f +  cos'P(5,0i 


(58) 


(59) 


The  global  x(s,t)  mdy(s,t)  are  related  to  the  slope  angle  ^(Sjt)  by 


^  =  cosmO 


(60) 


5 

x(5,  t)  =  x(0,  r)  +  J  cos^(  o;  t)d  <T 
0 

5 

y(5,O  =  >'(0,O  +  Jsin'P(ci;O^^<^  (61) 

0 

If  the  second  time  derivatives  of  x(s,t)  and^^s,^)  in  Eqs.  (60)  and  (61)  are  taken,  we  see  the 
transcendental  integration  of  the  unknown  angle  ¥(s,t)  over  space  will  include  terms  with 
the  first  two  time  derivatives  of  ^(s.t).  This  unknown  angle  ^(s,!)  is  also  contained  within 
trigonometric  fiinctions.  The  resulting  x  andj^  accelerations  are  dotted  with  the  unit  nprmal 
and  then  integrated  over  space  s  again  to  form  the  kinetic  energy.  This  leads  to  a  very 
complicated  and  non-standard  integro-differential  equation  form  for  the  equation  of 
motion.  The  double  spatial  integrals  and  trigonometric  functions  in  the  kinetic  energy 
greatly  complicates  any  approximation  of 'P(s,t)  and  the  resulting  solution  process. 

4. 1 .2  Finite  Element  Approach 

Principle  of  Virtual  Work 

A  method  due  to  Epstein  Murray^  is  paraphrased  in  this  section.  They  developed  a 
novel  nonlinear  formulation  which  is  applied  to  problems  involving  large  deformations  of 
elastic  beams.  Their  formulation  is  developed  using  the  principle  of  virtual  work  and 
implemented  using  a  firiite  element  discretization  of  the  beam  displacements.  The 
deformation  of  a  beam  is  defined  using  two  variables,  axial  displacement  u,  and  transverse 
displacement  w.  Normal  sections  of  the  beam  are  assumed  to  remain  plane  undistorted  and 
normal  to  the  beam  axis  after  deformation.  In  the  equations  that  follow  the  symbol 
O' denotes  an  s  derivative.  The  position  vector  to  a  point  on  the  axis  of  the  beam  is 

illustrated  in  Figure  12. 


=  sin'F(5,r) 
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Figure  12.  Finite  element  variables 

The  position  vector  to  a  point  on  the  axis  of  the  beam  after  deformation  is  given  by 

^  =  (s-¥u)i +yvj  (62) 

where  5  is  the  position  to  the  point  before  deformation  (arclength  measured  along  the 
undeformed  beam).  In  Cartesian  coordinates  the  position  of  the  deformed  point  is, 
X  =  (s+u),  andy  =  w.  The  position  of  a  general  point  not  originally  on  the  beam  axis  after 
deformation  is  given  by 


R  =  r+  ae^ 


(63) 


The  two  dimensional  Green's  strain  tensor  is  defined  using  derivatives  of  the  position 
vector  R.  The  only  nonzero  component  of  the  strain  tensor  is  derived  in  terms  of  the  axial 
strain  e  and  curvature  k  and  the  distance  off  the  neutral  axis  a 


s=  (1  - 


2 


(64) 


It  is  assumed  that  normal  sections  remain  plane  and  undistorted  after  the  deformation  of  the 
beam.  The  internal  virtual  work  (IVW)  is  defined  by  the  following  volume  integral 


W  =  JJ  (xSsdAdl 
I  A 


'(65) 


where  8  is  the  2-D  Green  Strain  measure,  a  is  stress  associated  with  8  and  /  is  the  length 
of  the  beam  axis  after  deformation.  The  variation  of  the  strain  component  8  is  taken  and 
substituted  into  the  above  equation.  The  integral  over  the  deformed  length  is  converted  to 
the  integral  over  the  original  length  using  the  definition  of  axial  strain.  After  further 
algebraic  manipulations,  the  terms  are  grouped  into  stress  resultants  and  the  following 
constitutive  equations  are  used 

Modified  internal  force  =  EAe  (66) 

Internal  moment  =  EI^  ~Z 

The  final  result  is 

WW  =  \  (EAe  5e  +  El  ^ 

0 


The  Green  Strain  measure  of  the  axis  is  given  by 


where  is  the  tangent  vector  along  the  beam  axis  after  deformation  and  i  is  the  base 
vector  along  the  undeformed  beam  axis.  The  tangent  vector  equals  the  derivative  of  the 
position  vector  given  in  Eq.  (62).  This  derivative  is  taken  and  the  axial  strain  e  is  written 
in  terms  of  the  displacements 


e  =  u'+|(»’2+w'2)  (70) 

The  Green  Strain  measure  is  also  used  to  develop  the  relationship  between  the  deformed 
length  and  the  undeformed  length 

dl  =  V2e  +  i  ds  (71) 


The  modified  curvature  x  is  developed  using  the  derivative  of  the  normal  vector  e„  and 
can  be  written  in  terms  of  the  displacements  u,  w  and  their  derivatives 

^  =  >i;"(l  +  m')  -  w'u"  (72) 


It  is  assumed  that  either  end  of  the  beam  can  be  subjected  to  an  external  force  T  and  S  and 
a  moment  M.  The  external  virtual  work  (EVW)  expression  for  these  boundary  conditions 
would  be 


EVWends  =  [tSu  +  SS^>+M^) 


s=I 

s=0 


(73) 
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•  where  the  angle  'F  =  tan  ’7-^  and  its  variation  is  ^  M\  +  u')S^'-w' Su'\ 

This  is  the  same  angle  as  the  slope  angle  'P  discussed  in  the  static  analysis.  For  small 
strains  the  denominator  (2e+l)  is  taken  equal  to  1. 

Finite  Element  Discretization 

A  set  of  algebraic  equations  for  Eq.  (68)  and  Eq.  (73)  is  developed  using  the  finite 
element  method  of  dividing  the  beam  into  elements  and  describing  the  displacement 

#  coordinates  u  and  w  in  terms  of  shape  functions  and  nodal  (endpoint)  displacements.  For 
each  element  the  displacements  are 


u(s,t)  =  ^fi(s)ui(t) 


i=l 


w(s,t)  =  'Yafi(s)wi(t) 
i=l 


(74) 

(75) 


The  shape  functions  chosen  are  the  cubic  polynomial  functions. 


hj 


f2=h 


h. 

,\2 


-i- 

<hj 


V 


h. 


\hj 


J 


(76a) 

(76b) 

(76c) 

(76d) 
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where  h  is  the  length  of  the  Ith  element.  The  beam  elements,  nodal  displacements  and  the 
numbering  scheme  for  u  and  w  are  illustrated  in  Figure  13. 
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11 
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W2=W6 
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Figure  13.  Element,  node  and  displacement  numbering 


The  displacement  discretizations  are  used  in  the  equations  for  the  strain  measure  e,  Eq.  (70) 
and  the  curvature  Eq.  (72)  and  their  first  variations.  This  is  substituted  into  the 
expressions  for  IVW,  Eq.  (68)  and  EVW,  Eq.  (73).  The  resulting  equations  are 

IVW  =  s  A /  {N*j  SU2^i-\)+j  +  N**k  (^2(/-l)+ j  ^2il-l)+k  +  ^2(/-1)+7  ^2(7-1)+^ )+ 

Mj  SW2(^i-i)+ j  +  {y2(i-i)+k  <^i^2(/-l)+ J  +  ^U2(I-l)+k)^  {^Jk  ~  ^Ig  )} 


where 

U2(,I-\)+kPjk  +2  (^2(J-l)+)t^2(/-l)+m  +  ^2(7-1)+* ^2(/-1)+/m)^7^ 
=E^2{I-V)+mPjkm  +^(^2(7-l)+m^^2(7-l)+K  +^2(7-1)+ot^2(7-1)+«)^7^« 
M*j  =  £7j^ir2(7-l)+)t-^;it  +^2(7-l)+ifc^2(7-l)+»i-0*»»«  ~^2iI-\)+k^2{I-l)+TnEjmk\ 

Mjl  =  +^(7-l)+?K^2(7-l)+n^ywfoj  ”^^(7-l)+m^2(7-l)+K^;«OTA:_ 
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The  upper  case  subscript  I  represents  the  element  number.  The  lower  case  subscripts  y,  k, 
m,  and  n  appearing  as  indices  in  two  adjacent  terms  represent  summations  for  each  index 

from  1-4.  For  example  the  term  represents  the  summation  of 

16  terms  with  m=l,4  and  n=l,4.  The  F  terms  represent  integral  evaluations  of  the  shape 
fiinctions  first  and  second  derivatives.  For  example 


where  c, 


s 

H' 


The  external  virtual  work  expression  is  brought  to  the  algebraic  form 


EVW  = 


^iXj SU2j+i+Yj SW2J+1  +  Mj^  + 1/21+2)^21 +2  ^7+2^27+2]} 

7=0^ 


where  7/  and  X/  represent  the  transverse  force  and  longitudinal  force  respectively  at  node 
I,  and  Mi  represents  the  moment  at  node  /.  The  equilibrium  equations  are  determined  by 
assuming  all  but  one  variation  of  the  nodal  coordinates  vanish.  There  are  two  contributions 
from  adjacent  elements  for  each  nodal  equilibrium  equation.  There  are  four  coordinates  at 
each  node  so  four  equations  are  generated  at  each  node.  The  following  four  equations 
apply  at  each  node 


8W2„.+i: 


h 


m  +  l 


He  HeHe  f 

M 3  +  Nj^  - 1)  +  7  +  -  ^j3  p2{m  -l)  +  j 


=  7 


m 


(79) 


5W2m+2‘ 


[He  He*  {  ♦  *  **^ 

M2  +  ^2m  +  j  ~  ^ j2  )^2m  +  y  J  + 

[H«  **  if^He*  **^  ,,  Arr  A 

Af4  +  AT ^-4 rr2(„ _  1) y  +  [ Af4y  - Mj^  p2(m-\)*  j\  =  ^m^*^2m-^2) 


(80) 


5U2m+l  : 


h 


m  +  \ 


He  **  1^  HeHe  ,  ,HeHe^^^ 

^1  ^J1  ^2m  +  7  v^yi  ”  7  J  2m  +  j 


h 


m 


He  HeHe 


He* 


**>| 


^3  +  ^ys ^2{m  - 1)  +  y  l^y3  “  ^3y  f2{m  -l)  +  j 


(81) 


=  X, 


m 


5U2m+2  : 


[He  He  * 

A^2  +  iN^y2  ^2/m  +  j  +  l,^y2  "  ^2j  J^2m  +  j_ 

[*  *  * 

N^+Nj^U2{jn-\)  + 


{  ** 

**" 

\  1 

[Mj,  - 

- 1 

+ 

I 

s 

(82) 


=  -^m^2m  +  2 


The  subscript  m  in  the  above  expressions  represents  a  node  number  from  0  to  N.  The 
equations  generated  are  nonlinear  due  to  the  multiplication  of  the  U  and  W  variables. 
Therefore  a  conventional  linear  stiffness  matrix  cannot  be  formed.  The  equations  for  all 
the  nodes  can  be  grouped  in  vector  form 


{NLTiq)}  =  {Q{q)} 

where  =  ^\iW2,U\,U2,W2,,W^,U2,,U . ,I^Ar+ljl^A/^+2»^2Ar+l»^2Ar+2} 


(83) 


The  virtual  work  done  by  the  internal  forces  are  nonlinear  terms  and  are  denoted  simply  by 
the  NLT  vector.  The  Q  vector  represents  the  nonconservative  work  vector.  The  only 
solution  approach  known  for  integrating  this  nonlinear  system  is  direct  numerical 
integration.  We  adopt  a  4*  order  Runge  Kutta  integration  scheme  to  solve  for  the  nodal 
forces  and  moments  given  a  shape  q.  Since  the  usual  problem  is  to  find  the  static  shape 


given  a  force  distribution,  a  Newton’s  method  is  used  to  solve  for  the  roots  of  the  set  of 
nonlinear  equations. 


Let  b{q)  =  {NLT{q)]-{Q{q)} 


then  a  Taylor  series  expansion  about  the  nth  iteration  q^  is 


^i(?)k+1  +■ 

=hi^)n  + 


^4N  . 


[<l4Nn+\-^^Nn 
(?4Ar„+i  -^4iV 


) 


(84) 


^47v(9)k+1  =*4iv(^)n  + 


‘^‘^4iV 

^4N 


[q4Nn+\--^4Nr) 


(85) 


This  requires  first  order  derivatives  of  each  nonlinear  equation  with  respect  to  each  variable 
in  q.  The  derivatives  of  the  equations  with  respect  to  these  variables  can  be  grouped  to 
form  the  Jacobian  matrix.  The  Jacobian  matrix  is  straightforward  to  program  but  is  a 
computationally  time  consuming  task  due  to  all  the  summations  in  the  multiplication  of  the 
coordinates. 


Hamilton ’s  Formulation 

The  equations  of  motion  are  formulated  using  Hamilton’s  principle  for  ease  of 
introducing  the  kinetic  energy  terms  and  allowing  the  conservation  of  energy  to  be  verified. 
Hamilton’s  principle  divides  the  energy  into  a  potential  energy  V,  a  kinetic  energy  T,  and  a 
nonconservative  work  term  W„c. 


The  potential  energy  function  due  to  the  strain  energy  of  axial  deformation  and 
transverse  deflection  is 


The  variation  of  the  potential  function  is  consistent  with  the  variation  of  internal  virtual 
work  developed  in  Eq.  (77).  The  resulting  nodal  equations  Eqs.(79)  -  (82)  can  be  directly 
imported  to  this  development. 

The  kinetic  energy  expression  can  be  derived  using  the  position  vector 
r  =  (s  +  u)i_  +  wj_ 

Since  the  components  are  in  inertial  coordinates,  the  velocity  is  simply 


r  =  ui_  +  wj 

The  kinetic  energy  expression  is 

T  =  (90) 

where  is  the  mass/length  of  the  beam.  The  variation  of  the  kinetic  energy  is 


54 


I 

5  J’=-Jp  [m5m  +  w6w]<3?5  (91) 

0 

Discretizing  the  variables  u  and  w  using  the  same  finite  element  shape  functions  in  Eqs. 
(76a)  -  (76d)  yields 


N 

7=1 


^^2(7-1)+ jP ik^2iI-\)+k  +  Pjk  <^2{I-l)+k 


(92) 


Again  the  lower  case  subscripts  represent  summations  from  1-4  so  the  group  expressions 
for  U  and  W  represents  the  sum  of  16  terms.  The  Fjk  term  represents  the  integral 
evaluations  of  the  shape  functions  listed  in  Eq.  (76) 

Fjk=\fjfk<^^  ^  =  J 

The  dynamic  equations  of  motion  are  devloped  from  the  kinetic  energy  variations  assuming 
all  but  one  of  the  node  variations  vanish.  There  are  four  equations  generated  for  each  node 
and  in  each  equation  there  is  a  contribution  term  from  each  element  adjoining  the  node. 


6W2m+i  : 

(94) 

5U2m+l  : 

K+lp[^2m+ jpjl]+  hm/^2{m-l)+jFj3\ 

(95) 

5W2m+2  : 

j^j2  j  /^[ 

+7-^74] 

(96) 

8U2m+2  : 

^ffj+1  /^[^2ff2+ j^J2 

y2(7«-l)+7^74] 

(97) 

The  lower  case  subscript  j  represents  a  summation  index  from  1-4  and  m  represents  a  node 
number  from  0  to  N,  An  element  mass  matrix  can  be  constructed  by  writing  these 
expressions  in  matrix  form  and  a  global  mass  matrix  is  developed  by  combining  the 
element  matrices.  Any  point  masses  can  be  added  to  the  appropriate  location  in  the  mass 
matrix.  Combining  these  results  with  the  potential  energy  results,  the  dynamic  equations  of 
motion  are 

[M]£  +  {iVLn^)}  =  [QiqS] 


For  this  nonlinear  system,  the  only  solution  approach  known  is  some  form  of  direct 
numerical  integration.  Numerical  integration  is  most  conveniently  done  in  terms  of  first 
order  equations,  which  requires  Eq.  (98)  be  rewntten  in  state  space  form. 


then 


\M]-^{Q{q,q)-NLm] 


(99) 


A  4*  order  Runge  Kutta  integration  scheme  is  used  to  solve  for  the  shape  and  velocity  over 
time. 


Friction  Forces 

The  nonconservative  forces  and  moments  are  included  in  the  generalized  force 
vector  Q.  Only  fiiction  forces  are  included  as  non-conservative  forces  in  this  study.  The 
generalized  force  is  developed  for  the  sliding  fiiction  between  a  mass  at  a  point  on  the 
beam  system  and  the  supporting  surface.  The  variation  of  the  nonconservative  work  is 


and  the  friction  force  is 

Fj-j.=-jUNsign{tpf)  where  N  =  mass  •  gravity  (101) 

The  variation  and  time  derivative  of  r  is  taken  from  Eq.  (89).  The  nonconservative  work 
becomes 

SWnc  =  -/^Nsign{upt"jSupt  -  jUN sign{wpf)Swpt  (102) 

and  the  coefficients  of  the  variations  are  added  to  the  generalized  force  vector  Q  at  the 
appropriate  locations. 


Constrained  Boundary  Conditions 

Constraint  equations  are  introduced  to  impose  boundary  conditions  and  to  allow  the 
coordinates  associated  with  a  fixed  support  to  be  removed  from  the  solution  process.  For  a 
clamped  boundary  condition,  the  w,  w  ’  and  u  coordinates  are  constrained  to  be  zero.  The 
constraints  are  written  in  matrix  form 

'l  0  0  0  0  0  .  0 

0  1  0  0  0  0  .  0 

0  0  1  0  0  0  .  0 

Let  the  matrix  on  the  left  hand  side  be  represented  be  C.  Lagrange  multipliers  are 
introduced  and  used  to  multiply  the  constraint  variations.  This  allows  the  constraints  to  be 
included  in  the  equations  of  motion  as 


0 


(103) 
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[M]£  +  {iVL7’(£)}=  {0(£4)}+  (104) 

The  Lagrange  multipliers  X  are  determined  by  differentiating  the  constraint  equation, 
substituting  in  the  equation  of  motion  for  q  and  solving  for  X 


yZ=[cM"^c]  ^CM-^\NLT(q)-Q(q,q)} 


(105) 


The  expression  for  X  is  now  substituted  into  Eq.  (104)  and  the  revised  equations  of  motion 
become 

[M]  q  +  {iVLr(£)}  =  {2(£,£)}+  C^[CM“^C]  ^ CM"^ {NLTiq)  -  e(£,£)}  (106) 


The  state  space  form  of  this  differential  equation  becomes 


and 


^2 


[M] 


{e(?.£)  -  iVir(£)}{l  -  C^{CM-'C^)~' CM  ‘ 


The  Runge  Kutta  4th  order  integration  technique  is  used  to  solve  these  equations  for  the 
displacement  shape  over  time.  The  computer  code  developed  for  the  forward  solution  of 
the  nonlinear  finite  element  model  is  listed  in  reference  14. 

4.1.3  Model  Validation 

Several  methods  are  used  to  validate  the  nonlinear  finite  element  model.  These  include 
matching  the  static  deformation  shapes  obtained  from  the  arc  length  model,  matching 


results  from  classical  linear  models  with  small  deformations,  and  checking  the 
conservation  of  energy. 

Verification  Using  the  Arc  Length  Approach 

The  equilibrium  shapes  obtained  from  the  finite  element  approach  and  the  arc 
length  approach  are  compared  for  two  loading  conditions.  The  first  comparison  is  for  the 
cantilever  beam  system  acted  on  by  a  moment  at  the  free  end.  This  example  was  used  in 
Figure  6  for  comparing  the  arc  length  approach  to  the  exact  solution.  Figure  14  shows  the 
finite  element  solution  using  four  elements  and  the  arc  length  solution  for  various 
moments. 


Figure  14.  Comparison  of  deflections  due  to  tip  moment 
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The  finite  element  approach  agrees  well  with  the  arc  length  solution  for  even  larger 
deformations  caused  by  a  pure  bending  moment.  For  a  static’s  problem  with  a  pure 
bending  moment,  the  two  methods  are  theoretically  equivalent  to  within  die  finite  element 
discretization  and  convergence  errors.  If  the  number  of  elements  is  increased  to  eight,  the 
finite  element  method  results  are  graphically  identical  to  the  arc  length  results  for  these 
bending  moments.  This  indicates  rapid  convergence  of  the  spatial  discretization 
approximation  to  the  arc  length  solution  of  Chapter  III. 

The  second  validation  comparison  makes  use  of  the  two  beam  PACE  system  with 
beam  properties  listed  in  Table  1.  This  system  is  subjected  to  an  external  force  acting  on 
the  free  end.  The  force  has  transverse  and  longitudinal  components.  The  prototype 
mechanism  designed  by  I.  Romero  allows  the  static  tip  force  to  be  directed  anywhere 
within  the  plane  of  the  beam.  The  finite  element  approach  is  applied  to  a  variety  of  loading 
configurations.  For  the  loading  condition  above  the  buckling  force,  there  are  two  solution 
shapes  obtained  from  the  finite  element  model  depending  on  the  initial  conditions  used  to 
start  the  nonlinear  iteration.  This  lack  of  uniqueness  occurs  only  above  the  critical  buckling 
load.  This  is  demonstrated  in  Figure  15.  The  lines  labeled  a  md  b  illustrate  the  two 
solutions  obtained  for  the  same  transverse  and  longitudinal  forces.  The  net  force  applied  on 
line  b  is  an  axial  compressive  type  force,  where  the  net  force  on  line  a  is  more  of  a 
transverse  compressive  force.  They  are  not  mirror  images  of  each  other  because  of  the 
effects  of  the  axial  versus  transverse  loading.  The  lines  a  and  a  do  represent  mirror 
images  where  the  longitudinal  force  is  the  same  and  the  transverse  forces  are  of  opposite 
sign.  Line  h  ’  is  the  other  nonlinear  solution  corresponding  to  the  same  loading  as  line  a 
It  is  a  mirror  image  of  line  b.  The  critical  buckling  load  is  calculated  to  be  12.7  N.  and  the 
loads  represented  in  Figure  15  have  a  longitudinal  force  of -15  N  and  a  transverse  force  of 


+.9  Nor -.9  N. 
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Figure  15.  Finite  element  solutions  for  tip  force 


^  It  is  difficult  to  compare  the  shapes  obtained  from  a  force  boundary  condition  on  the  finite 

element  model  to  the  arc  length  results  since  the  arc  length  modeling  approach  of 
Chapter  rH  does  not  include  the  effects  of  axial  deformation.  Therefore  the  models  are 
compared  in  Figure  16  for  small  transverse  loading  configurations.  It  is  expected  that  the 
•  two  nonlinear  solutions  will  converge  toward  each  other  and  the  classical  Euler  Bernoulli 

analytical  solution  for  small  deformations. 


1 


The  three  methods  agree  well  for  small  transverse  loading  but  as  expected,  the  axial  effects 
become  increasingly  significant  for  larger  deformations  and  the  methods  start  to  diverge. 

Static  Verification  for  Linear  Systems 

The  nonlinear  finite  element  model  is  checked  to  make  sure  it  recovers  the  classical 
results  of  a  beam  undergoing  axial  deformation.  The  equilibrium  equation  for  a  linearly 
elastic  bar  is 

— 1  =  0  (108) 

dx\  dx) 

If  the  bar  is  cantilevered  on  one  end  and  subjected  to  an  axial  load  p  on  the  other,  the 
analytical  solution  is 


The  axial  strain  du/dx  is  determined  and  compared  to  the  results  of  the  finite  element 
solution  for  various  forces.  The  results  are  listed  in  Table  6. 


Table  6.  Axial  strain  com] 

parison 

Axial  Force 

Strain  -  Finite  El 

Strain  -  Analytic 

-5 

-.364500e-06 

-.364500e-06 

-12 

-.874800e-06 

-.874799e-06 

-15 

-1.093501e-06 

-1.093499e-06 

-20 

-1.458002e-06 

-1.457999e-06 

The  axial  strain  produced  from  the  nonlinear  finite  element  model  matches  the  analytical 
solution  for  various  compressive  loads. 

The  nonlinear  finite  element  model  is  also  checked  against  the  analytical  solution 
for  a  beam  with  a  transverse  load.  The  equilibrium  equation  for  the  transverse  deflection 
of  an  Euler  Bernoulli  beam  is 


,2  r 


dx'‘ 


El 


dx^ 


=  0 


(110) 


The  analytical  solution  for  a  cantilever  beam  of  length  L  subjected  to  a  tip  force  F  is 
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#  y(x)^^(3L-x)  (111) 

The  two  beam  PACE  system  is  subjected  to  various  tip  forces  and  the  tip  deflections  are 

#  shown  in  Figure  16.  Both  the  finite  element  and  arc  length  methods  agree  with  the 
analytical  solution  for  small  loads  but  as  expected  do  not  agree  with  it  for  larger 
deformations  since  the  Euler  Bernoulli  beam  does  not  include  foreshortening  effects. 

^  Linear  Finite  Element  Dynamic  Solution 

A  linear  finite  element  model  is  constructed  for  the  PACE  two  beam  system.  The 
mathematical  model  is  developed  using  the  strain  energy  and  kinetic  energy  of  a>eam 
undergoing  axial  and  transverse  vibration*^.  The  potential  energy  function  includes  terms 

#  for  the  strain  energy  due  to  axial  motion  w  and  the  bending  energy  due  to  a  transverse 
deflection  u  of  an  Euler  Bernoulli  beam 


1 

V  =  -\  [EA{uf  +  El{w"f)ds 

The  kinetic  energy  due  to  axial  and  transverse  displacements  is 


r  =  —  j  pa{u^  + 
^0 


(112) 


(113) 


*  The  beam  is  divided  into  elements  and  the  displacement  coordinates  u  and  w  are 

approximated  using  the  element  end  displacements  and  shape  functions.  The  cubic 
polynomial  shape  functions  listed  in  Eqs.  (76a)  -  (76d)  and  the  four  endpoint  coordinates  in 
Eq.  (74)  are  used  for  the  transverse  displacement  w.  The  following  coordinate 
®  approximation  and  two  linear  shape  functions  are  used  for  the  axial  displacement  u 
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2 

i=l 

fl(s)  =  l-^  (114) 

The  element  mass  and  stiffness  matrices  are  constructed  using  the  element  approximations 
in  the  potential  and  kinetic  energy  expressions  and  evaluating  the  integrals  of  the  shape 
functions.  The  element  matrices  are  combined  to  give  the  global  mass  and  stiffness 
matrices  and  Lagrange’s  equations  are  used  to  derive  the  following  equation  of  motion 

Mq  +  Kq^QiqS 

The  analysis  of  this  linear  model  is  used  to  help  verify  the  nonlinear  model  and  used  to 
give  several  insights  to  the  expected  behavior.  The  natural  frequencies  and  mode  shapes 
are  determined  for  this  system  assuming  a  free  vibration  response.  For  the  PACE  two 
beam  system  with  eight  elements,  there  are  24  natural  frequencies.  Table  7  lists  some  of 
the  natural  frequencies  of  interest. 


The  lower  frequencies  correspond  to  transverse  vibration  and  the  first  three  have  periods  of 
2.6  seconds,  0.57  seconds  and  0.061  seconds.  The  higher  frequencies  correspond  to  axial 
vibrations  and  the  periods  range  from  3.87e-04  seconds  to  0.657e-04  seconds.  The  mode 
shapes  0=[(t)i,  <l)n]  and  frequencies  co,  are  used  to  transform  the  coordinates  q  to  modal 

coordinates  q  and  uncouple  the  linear  model  equations. 

Let  £(0  =  0  7(0  (116) 

then  the  analytical  solution  becomes 

9  (0  =  O_7(0)  cos(  <y0  (11^) 

The  equations  are  solved  analytically  for  a  free  vibration  problem  starting  from  an  initial 
deformation  shape  and  the  results  are  compared  to  the  nonlinear  finite  element  solution.  A 
numerical  solution  of  Eq.  (115)  is  also  determined  and  used  to  check  the  Runge  Kutta 
integration  algorithm.  The  shape  used  for  the  initial  deformation  corresponds  to  the  static 
transverse  load  of  F=8  shown  in  Figure  16  for  the  nonlinear  finite  element  solution  is. 
This  shape  is  used  instead  of  the  static  solution  to  Eq.  (115)  since  the  same  starting  shape  is 
needed  to  compare  the  free  vibration  response.  Table  8  compares  the  tip  transverse 
deflections  over  the  time  period  2.6  seconds  for  the  three  methods. 


Table  8.  Transverse  tip  deflections  at  given  time  intervals 

Time 

Analytical-Linear  FE 

Integration-Linear  FE 

Nonlinear  FE 

.1  sec 

.6037445  m 

.6037453  m 

.6046646  m 

.2  sec 

.5405400  m 

.5405403  m 

.5424933  m 

1.0  sec 

-.4671448  m 

-.4671460  m 

-.4792404  m 

2.0  sec 

.0799039  m 

.0799040  m 

.0957824  m 

2.6  sec 

.6083068  m 

.6083096  m 

.6071464  m 

The  numerical  solution  compares  closely  with  the  analytical  solution.  The  nonlinear  finite 
element  solution  agrees  with  both  linear  solutions  to  within  small  errors  in  the  third  digit. 
This  comparison  is  shown  graphically  in  Figure  17. 


Figure  17.  Free  vibration  comparison  of  methods 


The  nonlinear  finite  element  method  closely  matches  the  transverse  deflections  of  the 
linear  finite  element  model.  This  graph  does  not  reflect  the  axial  displacements  of  the  tip 
which  the  nonlinear  model  accounts  for.  The  linear  finite  element  model  does  not  have  the 
coupling  between  the  transverse  and  axial  motions.  It  is  of  significance,  however,  to  note 


that  accurate  numerical  solutions  require  Runge  Kutta  time  steps  which  are  a  fraction  of  the 
shortest  period  participating  in  the  response. 

Conservation  of  Energy 

A  check  of  the  conservation  of  energy  for  the  system  is  another  way  to  verify  that 
the  results  of  the  dynamic  model  are  logical.  The  sum  of  the  potential  and  kinetic  energy 
over  time  should  be  constant  in  the  absence  of  nonconservative  forces.  If  there  are 
nonconservative  forces,  the  change  in  the  energy  over  time  should  equal  the  work  done  by 
these  forces^®.  The  potential  energy  expression  for  the  nonlinear  model  is  listed  in  Eq. 
(87),  and  the  expressions  for  axial  strain  and  curvature  are  listed  in  Eq.  (70)  and  Eq.  (72). 
The  finite  element  discretizations  are  introduced  at  this  level  and  the  potential  energy  is 
ivritten  as 

f^2(/  - 1)  +7  f^2(/  - 1) + A:  + 

\P'jkmn  (t^2(/-l)+ y*^2(/-l)+A:^2(7-l)+»i^2(/-l)+« 

+  2I/2(/-l)+y^2(/-l)+A:^2(/-l)+m^2(/-l)+K 
+  J^2(/-l)+ j^2{I-\)+k^2(I-\)+m^2{I-\)+n  ) 

+  ^2(7-1)+ j^'jkm  (^2(7-l)+it^2(7-l)+ffi  +  ^2(7-l)+A:^2(7-l)+w) 

W2{I-\)^j^2{I-\)+kP"  jk'^^2{I-\)+j^2{I-\)+k^2{I-\)+mP"  jkm 
+  ^(7-l)+y^(7-l)+A:^2(7-l)+ffi^2(7-l)+n(-^  jkmn'^^  mnjk) 
-^i^2(I-i)+j^2(I-\)+k^2(I-l)+m^"kmJ 

_+  ^2(7-1)+ y^2(7-l)+A:^2(7-l)+»!^2(7-l)  +ri^  kmjn  )  J  J 

(118) 

where  the  expressions  for  Fj']^  are  combinations  of  shape  fimction  derivatives  evaluated 

over  the  integral  0  to  1,  and  the  lower  case  indices  represent  summations  from  1  to  4.  The 
determination  of  the  kinetic  energy  is  simpler  since  the  mass  matrix  is  already  determined. 


1  ^  I 
^7=1 


EA 


+  E 


(119) 


T  =  \fMq 


The  dissipative  work  done  by  nonconservative  forces  is 

t 

Wnc=\F-tdt  (120) 

0 


and  taking  the  time  derivative,  the  power  (work  rate)  equation  becomes 

W„c=F-t  (121) 

This  equation  can  be  numerically  integrated  over  time  along  with  the  displacement 
variables.  Friction  forces  are  included  in  the  finite  element  model  and  that  force  is 
described  in  Eq.  (101).  The  derivative  of  the  nonconservative  work  for  a  fiiction  force  at 
location  pt  is 


Wnc  =  -^J^{uptsigniUpt)  +  WptSign{Wpt)^ 


(122) 


During  the  dynamic  motion,  the  numerical  value  of  this  work  should  equal  the  change  in 
the  total  energy  of  the  system.  This  work  integral  computed  energy  change  is  used  to 
check  the  instantaneous  energy  in  all  of  the  simulation  results. 


4.1.4  Model  Refinement 

Time  Step  Analysis 

The  integration  time  step  chosen  for  the  Runge  Kutta  numerical  solution  is  based 
on  an  analysis  of  a  similar  linear  problem.  The  linear  problem  is  derived  by  taking  the 
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dynamic  nonlinear  model  in  Eq.  (98),  separating  out  a  linear  stiffness  matrix  operating  on 
the  displacement  variables  and  grouping  the  remaining  terms  in  a  nonlinear  vector. 

Mq+Kq  +  NLTiq)  =  Q{q,q)  (123) 

The  nonlinear  vector  NLT  is  neglected  in  the  linear  analysis.  The  resulting  linear  model  is 
qualitatively  similar  to  the  linear  model  used  in  the  dynamic  verification  Eq.  (115).  The 
difference  between  the  models  is  due  to  using  the  cubic  shape  functions  in  Eqs.  (76a  -  76d) 
for  both  variables  u  and  w  in  Eq.(123).  The  natural  frequencies  and  mode  shapes  are 
determined  using  the  following  open  loop  eigenvalue  problem 

KS  =  ^  (124) 

-i  ‘  -i 

The  mode  shapes  are  normalized  using  the  mass  matrix 


—i  — f 

(125) 

=  i  =  l,2,....« 

(126) 

The  frequencies  are  determined  for  an  8  element  linear  model  of  the  two  beam  PACE 
system.  The  lowest  frequency  is  2.42882  rad/s  and  the  highest  computed  frequency  is 
180075.56  rad/s  which  correspond  to  the  periods  2.587  seconds,  and  0.3489e-04  seconds 
respectively.  Although  the  higher  computed  frequencies  are  poorly  converged  in  the  finite 
element  analysis,  their  small  participation  in  the  solution  will  lead  to  numerical  stability 
problems  if  ignored  altogether.  A  time  step  of  .le-04,  which  is  one  third  of  the  shortest 
period,  allows  a  stable  integration  of  the  nonlinear  finite  element  model  to  proceed. 
Significantly  larger  time  steps  cause  the  motion  to  diverge  due  to  numerical  instability 
typical  of  stiff  differential  equations.  The  time  step  must  be  small  enough  to  track  the  high 


frequency  dynamics  for  stability  of  the  numerical  method.  When  accurately  integrated,  the 
higher  frequencies  and  mode  shapes  do  not  have  a  large  affect  on  the  dynamic  response  as 
are  seen  graphically  in  the  results.  The  lowest  frequency  corresponds  to  a  period  of  2.587 
seconds.  Unfortunately  for  a  stable  integration  over  the  largest  period  a  large  number  of 
integration  steps  are  required.  This  motivated  a  search  for  an  order  reduction  method  to 
eliminate  the  high  frequency  coordinates. 

Number  of  Elements 

The  total  number  of  elements  used  in  the  nonlinear  finite  element  method  is  based 
on  an  analysis  of  the  linear  model  analogous  to  the  preceding  section  and  based^  on  a 
convergence  study  of  the  nonlinear  initial  deformation.  The  natural  frequencies  for  the 
similar  linear  PACE  model  are  determined  for  three  cases  of  different  element  numbers 
and  are  grouped  according  to  transverse  w  modes  or  longitudinal  u  modes.  Table  9  and 
Table  10  lists  some  of  the  frequencies  that  can  be  compared  for  different  numbers  ,  of 
elements  used  in  the  linear  model. 


1  Table  9.  Natural  frequencies  for  w  modes 

4  elements 

8  elements 

16  elements 

2.42882 

2.42882 

2.42882 

11.0030 

11.0030 

11.0030 

103.234 

102.710 

102.675 

165.081 

163.049 

162.889 

419.539 

373.378 

371.689 

586.853 

472.479 

469.017 

Table  10.  Natural  frequencies  for  u  modes 

4  elements 

8  elements 

16  elements 

1711.09 

1705.61 

mi.%1 

4784.38 

AlZlSfl 

4706.50 

21120.4 

20952.8 

20854.1 

24195.5 

23665.4 

23440.9 

42434.9 

41691.0 

41465.0 

The  lower  frequencies  of  the  8  element  model  agree  closely  with  the  16  element 
frequencies  and  the  higher  frequencies  agree  within  1%  or  less.  The  frequencies  for  the  4 
element  model  are  not  converged  as  well  with  the  16  element  frequencies  and  the 
differences  range  from  0%  to  25%.  The  differences  between  die  frequencies  are  more 
pronounced  for  the  transverse  modes. 

The  effect  of  the  number  of  elements  on  the  initial  nonlinear  deformation  shape  is 
analyzed.  The  nonlinear  static  deformation  shape  is  solved  using  Eq.  (83)  for  two  different 
tip  loading  conditions.  Table  1 1  lists  the  tip  displacement  variables  for  different  numbers 
of  elements  due  to  a  moment  applied  at  the  tip. 


Ta 

ble  11.  Beam  d 

ieflection  for  tip  moment  =  18.( 

0 

#  Elements 

Wtip 

W’tip 

Utip 

U’tip 

4 

1.0772644 

.7176522 

-1.027003 

-1.696657 

8 

1.0777037 

.7018290 

-1.043144 

-1.712356 

16 

1.0776955 

.7017399 

-1.043231 

-1.712428 

32 

1.0776949 

.7017383 

-1.043232 

-1.712429 

The  8  element  model  agrees  well  with  the  16  element  model.  The  values  of  the 
displacements  for  the  8  element  model  agree  to  within  4  digits  of  the  higher  element 
models.  Another  deformation  shape  is  analyzed  for  models  with  different  numbers  of 
elements.  In  this  case  a  static  transverse  force  is  applied  to  the  beam  tip.  The  resulting  tip 
displacements  are  shown  in  Table  12. 


1  Table  12.  Beam  deflection  for  transverse  tip  force  =  8.0 

#  Elements 

Wtip 

w’tip 

Utip 

U’tip 

4 

.620976 

.610819 

-.167227 

-.208298 

8 

.627036 

.611709 

-.169576 

-.208926 

16 

.628348 

.612011 

-.170073 

-.209151 

32 

.628433 

.612040 

-.170106 

-.209172 

The  displacements  for  all  the  different  element  models  agree  to  within  2  digits.  The 
biggest  improvement  appears  to  be  between  increasing  from  4  elements  to  8  elements.  For 
the  8  element  model  the  values  of  the  w  and  u  coordinates  are  within  .0003  to.0014  of  the 
coordinate  values  for  the  32  element  model. 

Based  on  the  results  of  the  linear  frequencies  and  the  static  deformation  shapes,  the 
8  element  model  is  chosen  for  the  simulation  study  and  parameter  recovery  method. 

Modal  Truncation 

For  the  8  element  nonlinear  model,  a  small  time  step  is  needed  for  the  stability  of 
the  integration  method.  In  an  effort  to  allow  larger  time  steps,  an  attempt  is  made  to 
truncate  the  finite  element  model.  This  order  reduction  is  motivated  by  a  conventional 
method  used  for  linear  problems.  The  most  common  order  reduction  method  for  structural 
dynamics  is  the  modal  truncation  technique.  The  linear  finite  element  model  is 
transformed  to  modal  coordinates  and  the  high  frequency  low  amplitude  modes  which  do 


not  contribute  significantly  to  the  dynamics  are  partitioned  out.  Skelton’s^^  idea  of  a  modal 
cost  function  is  used  to  rank  the  modes’  relative  importance  to  the  dynamic  solution 
depending,  in  essence,  on  the  fractional  distribution  of  the  system  energy  among  the 
modes.  The  mode  shapes  and  frequencies  obtained  in  Eq.  (124)  are  used  for  a  coordinate 
transformation  similar  to  the  one  in  Eq.  (116).  The  modal  coordinate  nonlinear  equation 
assumes  the  transformed  structure  below 


K 7+  (^^[nLT[(^  7)}  =0^0 

where  M  =  M  ^  =  1  (127) 

and  K  =  (^^K<^  =  diagiiOi,ca2>--^i) 

The  first  order  form  of  Eq.  (127)  is  integrated  in  time  using  a  subset  of  the  modal 
coordinates.  The  solution  is  then  transformed  back  to  spatial  coordinates.  The  cost 
function  used  to  rank  the  modal  coordinates  is  the  contribution  of  each  modal  coordinate  to 
the  energy  of  the  linear  system.  The  nonlinear  model  is  used  to  solve  for  the  dynamic 
response  but  the  linear  form  of  the  energy  is  used  for  ranking  the  modal  contributions.  The 
kinetic  and  potential  energy  in  modal  coordinates  is 


E  =  M  '/;+ ^dt  (128) 

h 

The  transverse  modes  dominate  the  kinetic  energy  contributions  and  the  axial  modes 
dominate  the  potential  energy  contributions.  The  first  mode  #1  with  a  frequency  of 
2.4288  rad/s  is  the  largest  contributor  to  the  kinetic  energy  and  mode  #15  with  a  frequency 
of  4732.06  rad/s  is  the  largest  contributor  to  the  potential  energy  of  the  model.  Using 
Skelton’s  modal  cost  idea,  fifteen  modes  are  chosen  to  be  the  subset  of  coordinates  which 
represent  most  of  the  dynamic  motion.  These  modes  are  listed  below 
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^Tr  =  [<^1»  ^5’  Ao->  ^15>  ^19>  ^0»  ^1>  ^3»  ^4>  <^25’  ^28]  (^29) 

The  last  mode  included  in  the  subset  has  a  frequency  of  115665.15  rad/s  and  period  of 

•  .000054  sec.  Truncating  above  this  mode  would  only  allow  the  integration  time  step  to  be 
doubled  from  .00001  sec  to  .00002  sec.  The  solution  to  a  free  vibration  case  is  compared 
for  the  truncated  model  and  the  full  order  model  for  a  time  interval  of  .65  sec.  The  initial 
deformation  shape  is  due  to  the  transverse  tip  force  of  8  N.  A  comparison  of  the  beam 

*  deformation  shapes  for  the  truncated  model  and  the  full  order  model  is  shown  in  Figure  18. 
The  shapes  at  two  times  are  compared. 


Figure  18.  Beam  shapes  for  truncated  and  full  order  models 


The  displacements  for  the  truncated  model  are  found  to  be  significantly  different  than  the 
displacements  of  the  foil  order  model.  In  comparing  the  output,  the  differences  between 
the  models  are  contained  in  the  nonlinear  vector  term  NLT(^'x\).  The  small  differences  in 
the  modal  coordinates  are  summed  many  times  and  multiplied  by  large  constants,  EAi, 
EA2.  The  modal  truncation  is  tested  for  neglecting  only  one  mode,  #8,  which  is  the  lowest 
contributor  for  the  transverse  deflection.  The  displacements  obtained  after  one  quarter 
cycle  agree  to  within  6  digits  to  the  displacements  of  the  foil  order  model.  The 
displacements  at  node  2  and  node  3  are  slightly  different  (agree  to  4  digits)  than  the  foil 
order  model.  Neglecting  mode  #8  appears  to  have  a  very  small  affect  on  the  motionof  the 
beam  system  and  appears  to  affect  the  middle  of  beam  1.  Other  transverse  deflection 
modes  which  have  a  small  cost  contribution  can  be  neglected,  however  this  does  not  allow 
a  larger  time  step  since  these  frequencies  are  much  smaller  than  the  axial  mode 
frequencies.  Modal  truncation  is  tested  for  neglecting  the  three  highest  axial  modes,  #31, 
#32,  #33.  The  highest  frequency  remaining  in  the  model  is  142621.9  rad/s  and  the  time 
step  is  increased  to  .000015.  The  displacements  obtained  after  one  quarter  cycle  agree  only 
to  within  1  digit  to  the  displacements  of  the  foil  order  model.  The  increase  in  the  time  step 
does  not  justify  the  truncation  of  the  highest  modes  since  the  results  are  not  in  very  good 
agreement,  especially  since  our  analysis  indicates  the  foil  order  model  is  needed  as  a 
reference  to  confidently  truncate  any  modes.  Thus,  unfortunately,  modal  truncation  does 
not  appear  attractive  as  a  means  to  enhance  the  solution  efficiency  or  accuracy  for  this 
particular  nonlinear  finite  element  model. 

Rotational  Inertia 

The  contributions  to  the  kinetic  energy  of  the  rotational  motion  of  the  beams  and 
the  tip  and  elbow  masses  are  analyzed  to  see  if  the  rotational  inertia  should  be  included  in 
the  model.  The  rotational  kinetic  energy  for  a  body  is 
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where  T  =  [(l  +  u')w'  -  w'm']  and  /  is  the  mass  moment  of  inertia  about  the  center  of  mass. 
The  moment  of  inertia  about  the  z  axis  for  a  beam  of  differential  cross  section  dx  and 
thickness  a  is 

(131) 

For  the  PACE  system,  the  thickness  of  each  beam  is  .00315976  m.  The  moment  of  inertia 
for  the  differential  cross  section  is  computed  to  be  4.426E-07  kg-m.  This  value  is 'small 
and  the  angular  velocity  is  typically  1.85  s'K  Therefore  the  rotational  inertia  for  the  beams 
are  not  included.  The  inertias  for  the  elbow  and  tip  masses  are  122.98E-04  kg-m  and 
15.244E-04  kg-m^  respectively.  The  angular  velocity  is  substituted  into  Eq.  (130)  and  the 
first  variation  is  taken.  The  results  are  added  to  the  nonlinear  finite  element  equations  for 
du  ’  and  6w '  at  the  elbow  and  tip  node  locations.  The  5u  ’  and  5w  ’  equations  can  be  written 
in  matrix  form 

^  iX  +  u')^  -w'(1  +  m0 

[-w'(l  +  w')  w'^ 

The  first  nonlinear  matrix  part  of  the  above  equation  is  added  to  the  mass  matrix  and  the 
second  nonlinear  vector  is  added  to  the  NLT  vector.  The  free  vibration  of  the  finite 
element  model  is  integrated  over  half  a  cycle  using  the  initial  condition  due  to  a  tip 
transverse  force  of  8  N.  The  velocities  of  the  elbow  and  tip  angles  range  from  -.23  s’^  to  - 
1.85  s'^  over  half  the  cycle.  The  rotational  kinetic  energy  is  small  due  to  these  small 
velocities  and  the  small  moments  of  inertia  The  tip  transverse  deflections  are  compared  in 
Figure  19  for  the  model  including  rotary  inertia  and  without  it. 


w']  ^2u'(w'il  +  u')-w'u’y\ 


«'  J  ^  2  w'(w'(l  +  u')-  w'w')j 


(132) 


The  deflections  of  the  beam  model  including  rotational  inertia  effects  follows  almost 
exactly  the  deformation  history  of  the  beam  model  without  these  effects.  The  rotational 
motion  caused  by  a  free  vibration  of  the  PACE  two  beam  system  does  not  significantly 
affect  the  results.  Furthermore  including  the  rotational  kinetic  energy  complicates  the 
model  since  it  results  in  a  nonconstant  and  nonlinear  mass  matrix.  Therefore  the  rotational 
dynamics  of  the  tip  and  elbow  masses  are  not  included  in  the  model. 


4.1.5  Simulations 


A  free  vibration  analysis  is  presented  for  the  nonlinear  finite  element  model  of  the 
PACE  two  beam  system.  This  is  the  base  solution  for  the  parameter  recovery  method.  The 


deformation  shape  associated  with  a  tip  force  of  8  N  is  used  as  the  initial  condition.  Each 
beam  is  divided  into  four  elements,  giving  a  total  of  ei^t  elements.  The  equation  of 
motion,  Eq.  (106),  for  the  system  is  integrated  in  time  using  a  Runge  Kutta  scheme.  The 
system  is  studied  with  the  tip  mass  and  elbow  mass  values  listed  in  Table  1  and  includes 
sliding  friction  effects.  The  values  for  the  coefficient  of  kinetic  friction  between  metal  and 
metal  surfaces  range  from  0.1 1  to  0.45.^*  For  the  PACE  two  beam  system,  the  supporting 
surface  is  an  air  bearing  table  and  the  values  of  the  friction  coefficients  will  be  significantly 
less.  A  friction  coefficient  of  |ie=0.01  is  used  for  the  elbow  mass  and  p,=0.02  is  used  for 
the  tip  mass.  Figure  20  and  Figure  21  show  the  time  history  of  the  Wtip  and  uup  deflections 
over  a  cycle. 
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The  period  for  the  nonlinear  model  of  2.6  seconds  is  shown  in  the  transverse  tip  deflections 
and  corresponds  to  the  lowest  period  for  the  similar  linear  problem.  The  damping  shown 
in  the  tip  deflection  is  due  to  the  work  done  by  friction.  Figure  22  shows  the  whole 
deformation  shape  of  the  two  beam  system  for  one  quarter  of  a  cycle. 
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X  distance  (m) 

Figure  22.  Beam  profiles  for  1/4  period 


An  energy  balance  is  done  at  each  time  step  to  verify  the  change  in  the  system’s  potential 
and  kinetic  energy  equals  the  work  done  by  the  friction  forces.  The  initial  energy  of  the 
system  is  2.25257  N-m.  Table  13  summarizes  the  energy  check. 


Table  13.  Energy  balance  for  a  period 

Time 

A  Energy 

Work 

0.1 

-.00821430 

-.00821428 

0.5 

-.17588910 

-.17588908 

1.0 

-.44118282 

-.44118282 

1.5 

-.54355840 

-.54355839 

2.0 

-.76443096 

-.76443094 

2.6 

-.94819115 

-.94819112 

The  change  in  energy  differs  negligibly  from  the  work  done  by  the  friction  forces  and 
verifies  the  programming  of  the  friction  forces  is  consistent  in  the  model. 

The  finite  element  model  is  also  integrated  over  time  without  the  friction  forces  to 
ensure  die  energy  of  the  system  is  constant.  The  potential  energy  is  largest  at  the  highest 
and  lowest  transverse  deflection  shapes  while  the  kinetic  energy  is  largest  when  die  beam 
passes  through  the  axis  where  the  transverse  deflection  is  zero.  The  value  of  the  total 
energy  is  constant  and  agrees  to  within  9  digits  over  time  to  the  initial  value  of 
2.25257295  N-m.  These  tests  verify  the  nonlinear  finite  element  model  is  correct  from  a 
conservation  of  energy  perspective  and  provides  a  basis  for  confidence  that  the  formulation 
is  valid  and  correctly  programmed. 

4.2  Free  Vibration  Parameter  Estimation 

The  beam  parameters  affecting  the  free  response  characteristics  of  the  system  are 
recovered  using  a  set  of  measurements.  There  are  several  parameters  that  may  not  be  well 
known  such  as  the  friction  coefficients  or  the  mass  per  unit  length  of  each  beam. 


4.2.1  Least  Squares  Estimation 


The  least  squares  analysis  developed  for  the  inverse  static  analysis  is  used  to 
recover  the  best  fit  parameters  for  the  free  response  motion.  The  free  vibration  parameters 
are  represented  by  a  vector  p.  Measurements  at  different  time  intervals  are  represented  by 

the  vector  ^meas'  The  least  square  differential  correction  formula  for  updating  the 
parameters  is  repeated  below 


^  Corneas  ^ model  ^Pn)) 


A 

where  Q  =  — r- 


(133) 


The  sensitivity  matrix  will  obviously  change  since  the  nonlinear  model  and  the  parameters 
are  different  for  the  free  response  solution.  The  sensitivity  matrix  is  developed  using  the 
parametric  differentiation  technique  . 


The  first  order  form  of  the  nonlinear  finite  element  model  given  in  Eq.  (107)  is 


V  =  f{t,v,p) 


(134) 


The  time  derivative  of  the  sensitivity  matrix  becomes 


dt  ^p 


^  ^  <^(0  ^p 


(135) 


This  matrix  has  both  an  explicit  and  implicit  dependence  on  the  parameters  p  since  the 
displacement  and  velocity  variables  v  depend  on  p.  This  differential  equation  is  integrated 
simultaneously  with  the  solution  to  Eq.  (107).  The  entire  sensitivity  matrix  is  integrated 
but  only  a  portion  is  used  in  Eq.  (133)  for  the  variables  in  v  that  correspond  to  the  variables 
that  are  measured. 


The  derivatives  in  Eq.  (135)  depend  on  the  parameters  chosen  for  recovery  and  are 
described  in  the  following  section  for  friction  coefficients  and  mass  density  parameters. 
For  the  case  where  the  nonconservative  force  vector  Q  is  only  dependent  on  friction  forces, 
Eq.  (107)  becomes 


/  =  v  = 


{Q-NLTiq)} 


(136) 


The  implicit  derivative  matrix  df/dv  can  be  decomposed  into  four  sub-matrices.  The  first 
(1,1)  sub-matrix  is  a  zero  matrix  and  the  second  (1,2)  sub-matrix  is  the  identity  matrix.  The 
fourth  (2,2)  sub-matrix  is  also  a  zero  matrix  since  there  is  no  velocity  dependence  in  the  Q 
or  TVLT  vectors.  The  third  (2,1)  sub-matrix  consists  of  the  derivative  of  the  NLT  vector  with 
respect  to  each  variable  (the  Jacobian  matrix  discussed  earlier)  and  multiplied  by  [M]”* 

and  the  factor  CMULTI  where  CMt/IT/ =  -l}.  Calculation  of 

the  Jacobian  matrix  is  a  time  consuming  operation  due  to  all  the  nonlinear  summation 
terms  in  the  model.  For  the  parameter  recovery  program,  the  Jacobian  must  be  calculated 
at  each  time  step  and  the  whole  solution  repeated  each  time  the  parameters  are  updated. 
This  translates  into  a  large  amount  of  computer  time  depending  on  how  long  a  time  period 
is  covered  by  the  measurements.  The  implicit  derivative  matrix  is  summarized  below 


/ 


0 

I 

[M]~\CMULTI] 

<mLT 

0 

(137) 


The  explicit  derivative  matrix  df/dp  is  different,  depending  on  die  type  of  parameter 
involved  in  the  differentiation.  For  the  friction  coefficient  parameters  the  explicit 
derivative  matrix  is 


fL 


[M]-l 


m 

[CMULTI]\- 


(138) 


The  derivative  dQ/d^  for  friction  forces  is  a  vector  of  zeros  with  the  value  -  at 

the  appropriate  locations.  The  explicit  dependence  of  the  model  /  on  the  mass  density 
parameters  is  in  the  mass  matrix.  The  folloAving  identity  is  needed  in  the  derivative 

(139) 

After  several  algebraic  steps  the  explicit  derivative  matrix  for  the  mass  density  parameters 
becomes 


' 

{0}  1 
_  _ 

'  [M]~\CMULTI] 

[CMULTI][NLT{q)-Q}\^ 

(140) 


All  these  implicit  and  explicit  derivatives  are  combined  and  inserted  into  Eq.  (135). 

A  set  of  parameters  is  used  in  the  simultaneous  numerical  solutions  of  Eq.  (107) 
and  Eq.  (135).  The  parameters  are  then  updated  using  Eq.  (133).  The  numerical  solution 
and  updating  process  is  repeated  until  the  changes  in  the  parameters  are  small*'^. 

Measurement  Errors 

Physical  measurements  of  the  tip  locations  over  time  will  have  some  uncertainty  or 
error  variance  associated  with  them.  These  errors  are  included  in  the  parameter  updating 
process  for  the  free  response  motion  using  the  weighted  least  squares  approach  developed 
in  the  inverse  static  analysis.  The  weight  matrix  will  consist  of  the  inverse  error  variance 


associated  with  each  measurement.  The  formula  for  updating  the  parameters  is  repeated 
below 

=Pn  C'*!) 

The  sensitivity  matrix  for  the  nonlinear  finite  element  model  variables  is  developed  in  the 
preceding  section.  The  numerical  solutions  of  the  nonlinear  finite  element  model  and  the 
derivative  matrices  are  repeated  with  the  updated  parameters  until  the  changes  in  the 
parameters  are  small. 

A  priori  Estimates 

Some  of  the  parameters  may  be  known  within  a  certain  confidence  level  and  this 
can  be  reflected  in  a  priori  estimates  and  errors  of  these  estimates.  The  uncertainties 
associated  with  these  a  priori  estimates  are  used  to  limit  the  changes  to  the  parameters  so 
unreasonable  values  don’t  result.  These  uncertainties  are  incorporated  using  a  covariance 
matrix.  The  a  priori  estimates  can  be  treated  as  an  additional  measurement  equation. 

\Papriori  ~  Pn\~  +1  ~  Pn]'^  priori  (1^2) 

The  equation  for  the  a  priori  error  is  combined  with  the  error  expression  for  the 
measurements  in  Eq.  ( 1 4 1 )  yielding 


^tot 


^meas  '  ^ model 
Papriori  Pn 


a 


{Pn+\-Pn} 


(143) 


The  variance  matrix  for  the  a  priori  estimates  is  assumed  here  to  have  the  diagonal 
(uncorrelated)  structure. 
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(144) 


where  is  the  estimated  standard  deviation  of  the  a  priori  estimate  of  parameter  1, 

etc.  The  overall  weight  matrix  is  made  by  combining  the  variance  matrix  W  of  the 

measurements  and  the  variance  matrix  T  of  the  a  priori  estimates:  A  = 


~w 

o' 

r 

The  least  squares  function  of  the  error  is  minimized  giving  the  update  equation  for  the 
parameter  vector  p 


.  .  r  .  . . -1  .  T  A  ^meas  '  ^ model  (Pn  ) 
I  Papnori  Pn 


where  A : 


a 

I 


(145) 


The  covariance  matrix  for  the  estimated  parameters  p  including  a  priori  information  is 

£^=(a^Aa)''. 


A  set  of  parameters  is  used  in  the  numerical  solutions  of  Eq.  (107)  and  Eq.  (135). 
The  parameters  are  updated  using  Eq.  (145)  and  the  numerical  solution  is  repeated  until  the 
parameters  are  converged.  The  variance  of  each  parameter  is  taken  from  the  diagonal  of 
matrix  and  is  used  to  decide  how  well  the  parameter  is  determined. 

Parameter  Uncertainties 

The  estimated  parameters  recovered  in  the  static  analysis  can  be  brought  forward 
along  with  their  associated  uncertainties  and  combined  with  the  recovery  of  parameters 


affecting  the  free  vibration  response.  This  will  allow  adjustments  to  be  made  to  the 
statically  determined  parameters  if  they  are  not  determined  well.  If  the  previously 
determined  parameters  are  well  converged,  this  new  recovery  program  will  maintain  these 
values  while  allowing  other  parameters  to  be  changed.  The  parameters  and  uncertainties 
recovered  from  the  free  vibration  analysis  are  optimistic  since  the  values  of  the  static 
parameters  were  assumed  to  be  known.  The  free  vibration  parameters  can  be  further 
adjusted  using  this  combined  analysis.  This  will  be  achieved  by  using  an  approach  similar 
to  the  one  taken  in  the  parameter  recovery  method  with  a  priori  estimates.  The  parameters 
recovered  by  the  static  deformation  analysis  will  be  represented  by  the  vector 

=  {eIi,EI2)-  The  parameters  recovered  by  the  free  vibration  analysis  will  be 
represented  by  the  vector  p fr  =  Ph Pi)-  The  statically  determined  parameters 

can  be  treated  as  an  additional  measurement  p^  and  an  additional  error  equation  is 
developed. 

-Ps}  =  [I]{AP5  } + ^5  (146) 

This  equation  is  combined  with  the  error  expression  for  the  free  vibration  measurements 
(similar  to  Eq.  (38))  yielding 


^tot  ~ 


^meas  '  ^ model  ^P)  I 
Ps-Ps  \ 
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(147) 


The  sensitivity  matrix,  which  is  the  derivatives  of  the  measurement  variables  with 
respect  to  the  static  parameters  must  be  computed  for  the  nonlinear  finite  element  model 
and  included  in  the  solution  process. 


An  overall  weight  matrix  is  made  by  combining  the  variance  matrix  W  of  the  free  vibration 


measurements  and  the  covariance  matrix  of  the  static  parameters:  A  = 


w 

0 

0 

Es. 

The  weighted  least  squares  minimization  3delds  an  update  for  the  parameters 


=  (A^AA)"^A^Aj 
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where 


A  = 
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This  will  allow  further  updates  to  the  static  and  free  vibration  parameters  to  account  Tor  all 
the  parameter  uncertainties  in  a  free  vibration  experiment. 


4.2.2  Results  -  Free  Vibration  Parameters 

The  set  of  measurements  used  to  recover  beam  properties  are  simulated  using  the 
numerical  solution  of  the  nonlinear  finite  element  model  described  in  Eq.  (107).  The  free 
vibration  response  is  computed  for  a  set  of  ‘true’  parameters.  The  response  is  due  to  an 
initial  deformation  shape  caused  by  a  tip  transverse  force  of  8  N.  To  illustrate  robustness, 
we  consider  a  relatively  short  time  interval  transient  solution.  Of  course  longer  time 
intervals  would  yield  greater  precision  in  the  estimated  parameters.  The  tip  deflections 
over  one  period  of  the  dominant  linear  mode  are  shown  in  Figure  20  and  in  Figure  2 1 .  The 
tip  Wf  and  Uf  locations  are  chosen  as  the  measured  variables  at  selected  time  intervals  over 
the  period  of  the  beam  deflection. 

A  subset  of  the  dynamic  parameters  is  chosen  to  illustrate  the  recovery  method. 
This  subset  consists  of  the  fiiction  coefficients  of  the  elbow  and  tip  masses  and  the  mass 
densities  of  the  beams  Pfr  =  Ph  Pi)  •  *true’  mass  densities  are  listed  in 

Table  1  and  the  ‘true’  fiiction  coefficients  are  pe=0.01  and  p,=0.02.  The  calculations  of  the 
tip  location  Uf  and  w^  at  five  selected  time  intervals  are  input  to  the  parameter  recovery 


program  giving  a  total  of  10  measurements.  The  sensitivity  matrix  Q  is  a  4x10  matrix  as 
follows 


^9(h) 

^W9(?i) 

Pwgitx) 

^p\ 

^Pl 

(tj) 

^W9(?l) 

Pug  (tj) 

ppx 

Pp2 

^W9(?2) 

<^^9(^2) 

Pwg(t2) 

Ppx 

^Pl 

Pu^{t2) 

Pug{t2) 

Ppi 

^Pl 

Pug(ts) 

Pugit^) 

Ppx 

^Pl  . 

Estimates  for  the  mass  densities  of  the  beams  and  the  friction  coefficient  are  used  to 
start  the  program.  The  numerical  solutions  of  Eq.  (107)  and  Eq.  (135)  are  solved  for  the 
time  intervals  specified  and  the  parameters  are  updated  using  Eq.  (133).  This  is  repeated 
until  the  changes  in  the  parameters  are  small.  Table  14  summarizes  the  true  and  recovered 
values  of  the  parameters  assuming  there  are  no  errors  in  the  measurements. 


Table  1 

.4.  Parameter  recovery  for  perfect  measurements 

Melbow 

M-tip 

Pi 

P2 

True  values 

0.010 

0.02 

.5320 

.5300 

Starting  estimate 

0.015 

0.01 

.5603 

.5437 

Values  after  1  iter. 

0.010 

0.02 

.5297 

.5300 

Values  after  3  iter. 

0.010 

0.02 

.5320 

.5300 

The  program  rapidly  converges  to  the  true  parameter  values.  Other  numerical  experiments 
indicate  a  large  domain  of  attraction,  with  very  reliable  and  practical  convergence.  The 


parameters  associated  with  the  second  beam  system  are  recovered  more  rapidly  than  those 
in  the  first  beam  system,  due  to  the  larger  kinetic  energy  effects  at  the  tip  of  the  system. 

Measurement  Errors 

The  simulated  measurements  from  the  free  response  using  the  ‘true’  parameters 
have  errors  introduced  using  a  Gaussian  random  number  generator.  For  each  measurement, 
an  error  is  created  from  a  random  number  with  zero  mean  and  an  associated  standard 
deviation  a.  The  errors  are  added  to  the  ‘true’  measurements  which  are  then  used  to 
determine  the  best  fit  parameters  in  the  model. 

The  tip  locations  and  at  nine  selected  time  intervals  that  span  one  quarter  of 
the  deformation  cycle  are  used  for  a  total  of  18  measurements.  The  standard  deviations 
used  in  the  errors  of  these  measurements  range  from  0.0015  to  0.0006  for  the  w  and  u 
coordinates.  The  larger  the  standard  deviation  value  the  less  weight  is  given  to  that 
measurement  in  the  program  and  the  less  it  affects  the  recovered  parameters.  Table  15 
summarizes  the  results  of  the  recovered  values  of  the  parameters. 


Table  15.  Parameter  recovery  over  1/4  cycle  wi 

ith  measurement  errors 

M-elbow 

P-tip 

Pi 

P2  . 

True  values 

0.0100 

0.0200 

0.5320 

0.5300 

Initial  estimate 

0.0150 

0.0100 

0.5603 

0.5437 

Recovered  values 

0.0109 

0.0232 

0.1289 

0.4993 

Recovered  a 

0.0039 

0.0115 

0.4099 

0.1634 

The  recovered  values  of  the  friction  coefficients  are  seen  to  be  within  16%  of  the  true 
values  and  within  one  estimated  standard  deviation  of  the  parameters.  The  mass  densities 
are  significantly  different  than  the  true  values  and  have  correspondingly  larger  standard 
deviations  indicating  they  are  not  well  determined  by  these  measurements.  These  relatively 


larger  errors  indicate  poor  observability  and  suggest  a  longer  time  interval  experiment,  or  a 
new  experiment  in  which  these  parameters  play  a  more  important  role.  A  longer  time 
interval  of  measurements  is  used  to  better  estimate  the  mass  densities.  Measurements  of  the 
w,  and  Ut  variables  for  eighteen  time  intervals  oyer  one  half  of  the  deformation  cycle  are 
used  giving  a  total  of  36  measurements.  Again  the  standard  deviations  of  the 
measurements  range  from  .0015  to  .0006.  Table  16  summarizes  the  parameters  and  the 
associated  errors  recovered  for  this  case. 


Table  16.  Parameter  recovery 

over  1/2  cycle  with  measurement  errors  .. 

H'clbow 

M'tip 

Pi 

P2 

True  values 

0.0100 

0.0200 

0.5320 

0.5300 

Initial  estimate 

0.0150 

0.0100 

0.5603 

0.5437 

Recovered  values 

0.0084 

0.0220 

0.4677 

0.5360 

Recovered  a 

0.0020 

0.0026 

0.1199 

0.0843 

The  values  for  the  friction  coefficients  estimated  are  much  closer  to  the  true  values  than  the 
coefficients  recovered  for  the  quarter  cycle  measurements.  The  small  error  variance  of 
these  parameters  indicate  they  are  fairly  well  determined.  The  recovered  value  of  the  first 
beam’s  mass  density  is  worse  than  the  initial  estimate  and  this  estimated  value  is  outside 
the  commonly  accepted  density  range  of  standard  aluminum  beams'^.  The  estimated 
standard  deviation  of  this  parameter  is  consistently  large  compared  to  the  value  of  the 
parameter,  so  at  least  the  estimation  process  “knows  it  doesn’t  know”.  This  indicates  the 
recovery  of  the  beam  1  density  parameter  cannot  be  trusted.  The  recovered  value  of  the 
second  beam’s  mass  density  is  close  to  the  true  value  and  is  an  improvement  over  the 
estimated  value.  However,  it  still  has  an  associated  standard  deviation  that  is  large  which 
would  indicate  the  value  of  the  parameter  is  not  well  determined  suggesting  that  this 
convergence  is  fortuitous.  Therefore,  the  densities  of  the  beams  are  weakly  observable 


and  there  is  little  confidence  in  these  recovered  parameter  values.  On  the  other  hand,  the 
other  parameters  can  be  estimated  from  a  very  short  transient  response. 

The  mass  densities  of  the  beams  are  clearly  difficult  to  recover  in  the  presence  of 
measurement  errors.  This  is  because  they  do  not  have  a  large  effect  on  the  free  response 
behavior  of  the  system  in  this  particular  experiment.  Table  17  lists  the  tip  locations  solved 
at  various  times  using  the  recovered  parameter  values  listed  in  Table  16.  The  measured  tip 
locations  input  to  the  model  and  the  tip  locations  obtained  using  the  true  parameters  are 
included  in  this  table  for  comparison. 


Table  17.  Effects  of  parameter  values  on  tip  transverse  deflections 

Time 

Wt  -  measured 

Wt  -  true  parameters 

Wt  -  rec.  parameters 

0.1 

.604670 

.605715 

.605819 

0.5 

.256023 

.255165 

.255321 

1.0 

-.413409 

-.413560 

-.413387 

1.3 

-.550037 

-.550117 

-.550549 

The  tip  deflections  using  the  recovered  parameter  values  are  closer  to  the  deflections 
predicted  by  the  true  parameter  values  than  the  measured  values.  This  is  due  to  the  effect 
of  the  number  of  measiurements  on  the  relationship  between  the  variance  of  the  parameters 
and  the  variance  of  the  measurements.  This  improvement  is  seen  even  with  the  unlikely 
values  of  the  beam  mass  densities  used  as  part  of  the  parameters.  Again,  the  apparent 
paradox  simply  reflects  the  truth  that  the  beam  densities  have  such  a  small  effect  on  the 
free  response  of  the  system  over  this  relatively  short  transient  motion.  If  there  were  no 
external  forces,  we  can  easily  show  that  the  mass  densities  would  not  be  uniquely 
determined,  so  the  observability  is  also  linked  to  the  smallness  of  the  friction  force. 


A  priori  Estimates 


The  starting  estimates  of  the  parameters  can  be  known  within  a  certain  tolerance 
and  this  information  can  be  used  to  limit  the  changes  in  the  parameters.  As  was  pointed 
out  in  the  last  section,  the  value  of  the  first  beam  mass  density  recovered  is  below  the 
known  density  range  for  aluminum  beams  having  that  shape.  A  priori  information  can  be 
used  to  limit  the  estimates  of  low  observability  parameters  to  be  consistent  with  an  a  priori 
estimate’s  uncertainty. 

The  tip  locations  for  eleven  time  intervals  over  one  quarter  cycle  of  the  period  are 
input  for  a  total  of  22  measurements  to  the  parameter  recovery  program.  The  standard 
deviations  used  for  the  measurements  range  firom  0.0015  to  0.0006.  A  priori  estimates  of 
the  parameters  are  used  as  the  starting  values  in  the  program.  The  a  priori  variance  values 
of  the  beam  mass  densities  and  fiiction  parameters  will  have  negligible  impact  on 
observable  parameters,  but  will  serve  to  hold  poorly  estimated  parameters  to  the 
neighborhood  of  the  a  priori  estimates  (consistent  with  a  priori  estimates).  The  converged 
covariance  matrix  will  correctly  reflect  the  total  information  content  of  both  the  a  priori 
information  and  the  measurements. 

The  numerical  solution  of  Eq.  (107)  and  Eq.  (135)  are  computed  to  the  final  time 
period  and  Eq.  (145)  is  used  to  update  the  parameters.  This  recovery  process  is  repeated 
until  the  change  in  the  parameter  values  is  small.  Table  18  summarizes  the  parameters 
recovered  and  the  converged  parameters’  standard  deviations. 


Table  18.  Parameter  recovery  with  a  priori  information 

M^lbow 

Ptip 

Pi 

P2 

True  values 

0.0100 

0.0200 

0.5320 

0.5300 

a  priori  estimate 

0.0150 

0.0100 

0.5603 

0.5437 

a  priori  a 

0.0075 

0.0150 

0.0532 

0.0530 

Recovered  values 

0.0081 

0.0203 

0.5486 

0.5557 

Recovered  a 

0.0019 

0.0042 

0.0526 

0.0481 

Including  the  a  priori  variance  matrix  stabilized  the  poorly  observed  densities  and 
indirectly  the  tip  mass  friction  coefficient.  The  recovered  value  of  ptip  is  closer  to  the  true 
value  than  the  recovered  value  of  ptip  with  only  the  effect  of  measurement  errors  included. 
The  value  recovered  for  the  elbow  mass  friction  coefficient  was  not  significantly  affected. 
The  a  priori  variance  values  did  limit  the  change  in  the  beam  mass  densities  and  kept  the 
recovered  parameters  closer  to  the  a  priori  as  well  as  the  true  values.  The  standard 
deviations  of  the  recovered  mass  densities  are  smaller  compared  to  those  obtained  using 
only  measurement  errors  but  they  remain  comparable  to  their  a  priori  values.  This 
indicates  a  large  range  of  densities  can  be  obtained  and  still  be  consistent  with  the 
measured  motion. 

Since  the  mass  densities  are  difficult  parameters  to  estimate,  they  are  removed  from 
the  parameter  estimation  and  only  the  friction  coefficients  are  estimated  using  the  a  priori 
covariance  matrix.  The  standard  deviations  of  the  friction  coefficients  allows  large 
changes  of  the  a  priori  values.  Measurements  of  the  tip  u  and  w  locations  are  taken  for  six 
time  intervals  and  similar  measurement  errors  are  used  as  listed  in  the  last  section.  Table 
19  lists  the  recovered  friction  parameters  and  their  associated  errors. 


1  Table  19.  Friction  coefficient  recovery  wit 

;h  a  priori  information 

M-elbow 

Mtip 

True  values 

0.0100 

0.0200 

a  priori  estimate 

0.0150 

0.0100 

a  priori  a 

0.0080 

0.0120 

Recovered  values 

0.0093 

0.0202 

Recovered  a 

0.0042 

0.0036 

The  values  recovered  for  the  friction  coefficients  using  the  a  priori  variance  matrix  are 
close  to  the  true  values  and  well  within  one  standard  deviation  of  the  true  values.  The 
a  priori  standard  deviations  of  the  friction  coefficients  are  doubled  to  allow  greater 
changes  to  the  parameter  values  and  the  parameter  recovery  program  is  repeated.  These 
results  are  shown  in  Table  20. 


Table  20.  Friction  coefficient  recovery  with  larger  a  priori  variance 

M-elbow 

M-tip 

True  values 

0.0100 

0.0200 

a  priori  estimate 

0.0150 

0.0100 

a  priori  a 

0.0160 

0.0240 

Recovered  values 

0.0067 

0.0224 

Recovered  a 

0.0049 

0.0041 

The  recovered  friction  coefficient  values  are  not  as  close  to  the  true  values  as  those 
recovered  for  the  more  limiting  variance  matrix.  They  are  within  one  standard  deviation  of 
the  true  values.  A  more  reasonable  process  would  be  to  repeat  this  experiment  with  a 
distribution  of  estimates  with  the  a  priori  statistics.  This  would  allow  one  to  confirm  the 


truth  that  the  estimation  process  is  statistically  consistent.  However,  any  given  experiment 
may  look  contrived  because  it  is  merely  a  sample  from  a  distribution  of  possibilities. 

The  recovered  friction  coefficients  are  used  in  the  program  to  determine  the  free 
response  deformation  shapes  of  the  beam  system.  The  tip  transverse  deflections  are 
compared  to  the  deflections  obtained  using  the  true  coefficients  and  compared  to  the 
measured  data  used  for  input  to  the  parameter  recovery  program.  This  comparison  is 
summarized  in  Table  2 1 . 


1  Table  21.  Effects  of  a  priori  determined  friction  parameters 

Time 

Wt  -  measured 

Wt  -  tme  parameters 

Wt  -  rec.  parameters 

0.1 

.606013 

.605715 

.605723 

0.2 

.546752 

.546841 

.546847 

0.3 

.463057 

.464102 

.464035 

0.4 

.368208 

.367896 

.367637 

The  beam  deformation  shapes  using  the  recovered  parameters  agree  well  with  the 
deformation  shapes  using  the  true  values.  There  is  3  to  4  digit  accuracy  between  the  tip 
values.  The  friction  coefficients  are  fairly  well-determined  parameters  and  we  have  the 
means  to  accommodate  consistently  any  a  priori  information  available. 


CHAPTER  V 


NONLINEAR  FORCED  VIBRATION  ANALYSIS 


5.1  Development 

Including  motor  effects  in  the  modeling  of  the  PACE  two  beam  system  is  the  last 
step  in  the  sequential  system  identification  process.  Additions  to  the  nonlinear  finite 
element  model  are  made  to  reflect  the  dynamics  of  a  motor  attached  to  the  structure. 


5.1.1  Motor  Modeling 

The  motor  operation  and  characteristics  needed  for  this  computer  simulation  are 
taken  from  hardware  commonly  used  in  the  Dynamics  and  Controls  Laboratory  of  the 
Aerospace  Engineering  Department  at  Texas  A&M  University^.  A  reaction  wheel  motor  is 
placed  on  the  tip  of  the  structure  to  suppress  end  effector  vibration.  A  torque  is  applied  on 
the  tip  of  the  structure  by  accelerating  the  reaction  wheel  with  a  DC  motor.  The  torque  of 
the  motor  is  a  function  of  the  voltage  sent  to  the  motor  by  the  power  amplifier^  ^  Since  the 
power  amplifier  can  be  operated  in  current  or  voltage  mode,  we  assume  the  power 
amplifier  is  operated  in  the  current  mode  which  avoids  including  the  dynamics  of  the  motor 
in  the  simulation.  This  is  because  our  motor  has  a  near  linear  torque/current  relationship. 
The  power  amplifier  uses  an  internal  control  law  to  adjust  the  voltage  supplied  to  the  motor 
to  follow  the  desired  current.  The  time  constant  of  the  power  amplifier  is  sufficiently  short 
( =  10'^  s)  that  the  response  time  is  neglected  in  this  study.  A  Maxon  DC  motor  outfitted 
with  a  reaction  wheel  is  sent  a  current  which  is  used  to  apply  a  commanded  control  torque. 

T  =  Kti  (149) 

where  Kt  is  the  torque  constant  of  the  motor  and  equals  0.03891  N-m/amp.  The  linear 
torque/current  relationship  is  accurate  to  less  than  1%  over  a  bandwidth  from  0  to  30  Hz 
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assuming  saturation  conditions  are  not  encountered.  A  control  law  to  suppress  end 
vibration  is  chosen  and  is  a  linear  feedback  damper  of  the  form 

T  =  -g^tip  (150) 

where  g  is  the  gain  and  is  the  angular  velocity  of  the  beam  tip. 

For  the  nonlinear  finite  element  model,  the  torque  is  included  in  the  nonconservative  work 
term 


SWnc^T5^,ip 


(151) 


Using  the  coordinate  system  established  in  Figure  12,  the  variation  of  the  tip  angle  is 


^tip  = 


|l  +M^/p  j  ^tip  ^tip  ^^tip  j 


(152) 


and  the  velocity  is 


^tip  ^tip  \^tip  ^tip  ^tip  j 


(153) 


Combining  Eqs.  (150)  -  (153)  gives  the  following  expression  for  the  nonconservative  work 


^nc  ^tip 

+  +  ^tip  \^tip'^iip  ~  {^tip  )  ^tip  j 


^^tip 


tip 


(154) 


The  terms  associated  with  Sw^p  and  Su'^pare  included  in  the  Q  vector  and  the 
numerical  solution  of  Eq.  (107)  is  used  to  determine  the  system  displacements  over  time. 

5.1.2  Model  Validation 

An  energy  balance  is  done  at  each  time  step  to  ensure  the  change  in  the  energy  of 
the  system  equals  the  work  done  by  the  motor  torque  and  the  friction  forces.  The  work 
done  by  the  motor  torque  is 


Taking  the  time  derivative  gives 
jr=Ti'=-g^^ 


(155) 


(156) 


This  equation  is  combined  with  the  work  equation  for  friction  forces  Eq.  (122)  and  the 
result  is  numerically  integrated  over  time  along  with  the  displacement  variables.  At  each 
time  step,  the  numerical  value  of  this  work  should  equal  the  change  in  the  total  energy  of 
the  system.  This  check  indicates  the  motor  torque  has  been  correctly  included  in  the  finite 
element  model. 

5.1.3  Simulations 

The  nonlinear  finite  element  model  is  integrated  forward  in  time  including  the 
effects  of  a  tip  motor.  The  control  gain  g  for  the  feedback  control  law  must  be  chosen  and 
can  be  determined  by  several  methods.  For  example,  the  finite  element  model  can  be 
linearized  and  control  design  methods  e.g.  a  pole  placement,  quadratic  regulator,  etc.  can 
be  used,  or  the  gain  can  be  chosen  by  trial  and  error  to  give  die  desired  damping.  For  this 
study,  the  trial  and  error  method  is  used  since  it  is  simpler  and  the  main  focus  of  this 


research  is  the  recovery  of  model  parameters  and  not  on  control  design.  Several  gain 
values  are  used  in  the  nonlinear  finite  element  model  and  the  resulting  deformation  shapes 
compared.  Figure  23  and  Figure  24  show  the  tip  w  and  u  displacements  for  different  values 
of  g. 


Figure  23.  Effect  of  gain  values  on  tip  transverse  displacements 
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Figure  24.  Effect  of  gain  values  on  tip  longitudinal  displacements 


The  higher  control  gain  of  2.0  has  the  largest  reductions  in  the  tip  deflections  and  brings 
the  beam  to  rest  quickly.  However  the  torque  needed  to  suppress  the  vibrations  is  large  and 
beyond  the  operating  bounds  of  a  Maxon  motor.  The  smaller  control  gain  damps  out  the 
tip  motions  and  generates  torques  within  the  bounds  of  the  motor.  The  highest  torque  this 
gain  requires  is  .35  N-m  which  corresponds  to  49  oz-in.  A  plot  of  the  tip  displacements 
and  the  torque  history  for  the  control  gain  of  0.2  is  shown  in  Figure  25. 
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Figure  25.  Control  and  tip  displacements  for  control  gain  of  0.2 


An  energy  balance  is  performed  at  each  integration  step  and  the  change  in  the  total  energy 
of  the  system  equals  the  change  in  the  work  done  by  the  motor  torque  and  the  friction 
forces.  These  values  match  to  eight  digits.  The  initial  energy  of  the  system  is 
2.2525729  N-m  and  after  one  period  the  energy  change  (and  work  done)  equals 
1.2967753  N-m. 

5.2  Forced  Vibration  Parameter  Estimation 

The  parameters  affecting  a  forced  response  are  related  to  the  motor  used  on  the 
system.  The  parameter  recovered  in  this  analysis  is  the  ratio  of  the  motor  torque  constant. 
If  the  torque  constant  is  not  well  known,  the  torque  generated  by  the  motor  is  incorrect  and 
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•  the  effects  on  the  recovered  deformation  shapes  are  used  to  recover  the  correct  torque 

constant,  .  If  the  torque  constant  is  incorrect,  ,  the  current  sent  to  the  motor  is 
incorrectly  computed  as 


(157) 


The  torque  applied  to  the  structure  is  not  the  desired  torque,  but  an  altered  value  related  to 
the  ratio  of  the  true  and  estimated  torque  constants. 


7’=r(-g'F)  where  r  =  -^  (158) 

This  equation  for  the  torque  is  used  in  the  development  of  the  nonconservative  work  done 
by  the  motor. 


(159) 


The  Sw'fip  and  terms  are  included  in  the  Q  vector  in  the  nonlinear  finite  element 
model  represented  by  Eq.  (107). 


5.2.1  Least  Squares  Estimation 

The  least  square  procedure  outlined  in  Chapter  IV  for  the  estimation  of  static 
•  parameters  is  used  to  estimate  the  motor  torque  constant  ratio.  The  true  torque  constant  is 

then  estimated  by  a  simple  multiplication  between  r  and  the  estimated  torque  constant. 
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To  test  these  ideas,  a  set  of  ‘true’  measurements  of  the  tip  Ut  and  Wt  deflections  is 
simulated  by  integrating  the  differential  equation  in  Eq.  (107)  including  the  motor  model 
from  Eq.  (159)  and  a  known  value  of  the  torque  constant  ratio.  The  recovery  of  the 
parameter  r  is  started  by  estimating  a  value  of  r  and  integrating  in  time  to  the  last 
measurement  interval.  The  parameter  value  is  updated  using  the  sensitivity  matrix  and  the 
difference  between  the  measured  and  model  tip  deflections  as  developed  previously. 

^  Corneas  model  (1^0) 

The  only  real  change  in  this  parameter  estimation  process  is  in  the  calculation  of  the 
sensitivity  matrix  Q .  In  this  case  the  sensitivity  matrix  is  composed  of  one  colunrn  since 
only  one  parameter  is  being  recovered  in  this  discussion.  The  rows  are  derivatives  of  the 
variables  unp  and  Wup  with  respect  to  r  and  evaluated  at  the  measurement  times.  The 
sensitivity  matrix  is  integrated  in  time  to  get  the  derivatives  at  the  measurement  times.  The 
general  form  of  the  differential  equation  for  the  sensitivity  matrix  is  the  same  as  in 
Eq.  (135).  However,  the  implicit  derivatives  are  different  since  including  motor  effects 
causes  a  dependence  on  the  tip  slopes  and  velocities  in  the  Q  vector.  The  implicit 
derivative  matrix  becomes 


0 

I 
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[MY\CMULTI] 
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(161) 


Since  Q  only  depends  on  w'^p ,  w/jp ,  w'tip ,  and  u'tip ,  the  derivative  matrices  ^Q/  are 
sparse  matrices  with  few  entrees.  The  explicit  derivative  matrix  is  also  different 

for  the  parameter  r.  It  becomes 
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These  derivative  matrices  are  used  in  Eq.  (135)  for  integrating  the  sensitivity  matrix  over 
time.  The  sensitivity  matrix  and  the  derivatives  included  in  this  matrix  are  checked  via 
finite  difference  calculations  at  each  time  step  and  there  is  good  agreement  between  them. 
This  verifies  the  correct  development  and  placement  of  the  derivatives  used  in  the 
integration  of  the  sensitivity  matrix.  An  estimate  of  the  motor  torque  constant  is  used  in 
the  simultaneous  numerical  solutions  of  Eq.  (107)  and  Eq.  (135)  with  Eq.  (159)  used  for 
the  motor  model.  The  parameters  are  then  updated  using  Eq.  (160).  The  numerical 
solution  and  updating  process  is  repeated  until  the  changes  in  the  parameters  are  small. 

Measurement  Errors 

The  parameter  recovery  of  the  motor  torque  constant  ratio  can  be  altered  to  include 
the  uncertainties  in  the  measurement  values.  The  same  process  is  followed  as  the  one  used 
in  the  recovery  of  the  static  parameters.  A  weight  matrix  is  included  in  the  least  square 
minimization  function  and  the  weights  used  are  the  reciprocal  of  the  error  variances  for 
each  measurement.  The  same  parameter  update  formula  is  developed  but  the  different 
sensitivity  matrix  is  used  in  the  equation  as  discussed  in  the  previous  section. 

'•«+!  ='•«  +(n^W5)->n^(r(Y„,„-Y„„*,(r„))  (163) 

The  measurement  errors  are  mapped  into  the  error  in  the  recovered  motor  parameter  using 
Ej.  =  ^  =  cr^ .  The  error  variance  of  the  motor  torque  constant  ratio  reflects  on 

how  well  determined  this  parameter  is. 


Parameter  Uncertainties 


The  estimated  parameters  and  uncertainties  from  the  previous  static  and  free 
vibration  analysis  can  be  carried  forward  and  combined  with  the  recovery  of  parameters 
affecting  the  forced  response.  Measurement  information  is  brought  forward  for  a  final  level 
of  parameter  updating.  This  new  recovery  program  will  take  advantage  of  well  converged 
parameters  and  allow  further  adjustments  to  be  made  for  ill  converged  parameters.  The 
parameters  are  grouped  by  the  recovery  method  used.  The  statically  determined  parameters 

are  Pg={EIi,EI'2),  the  parameters  recovered  by  the  free  vibration  analysis  are 
Pfr  -  P\>  stid  the  parameters  recovered  by  the  forced  response  are  Pf  =  (r) 

The  previously  determined  parameters  are  treated  as  additional  measurements 
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This  equation  is  combined  with  the  forced  vibration  measurement  error  expression  yielding 
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An  overall  weight  matrix  is  made  by  combining  the  variance  matrix  ^  of  the  forced 
vibration  measurements  and  the  covariance  matrix  of  the  static  parameters  and  the 
covariance  matrix  E^  of  the  free  vibration  parameters: 
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The  weighted  least  squares  minimization  yields  an  update  for  the  parameters 
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This  combined  parameter  update  allows  the  measurements  and  covariance  estimates  for  all 
the  parameters  to  influence  the  final  determination  of  the  static,  fi'ee  vibration,  and  forced 
vibration  model  parameters. 

5.2.2  Results  -  Forced  Vibration  Parameters 

A  set  of  ‘true’  measurements  of  the  tip  locations  are  simulated  by  integrating  the 
nonlinear  finite  element  model  Eq.  (107)  using  a  known  value  of  the  torque  constant  ratio 
in  Eq.  (159).  Any  known  value  of  the  torque  constant  ratio  can  be  used  and  a  value  of 
r  =  1  is  chosen  since  simulated  measurements  are  available  from  previous  work.  The 
values  of  w,  and  Ut  are  taken  at  five  time  intervals  over  one  quarter  cycle  for  a  total  of  ten 
measurements.  The  starting  estimate  of  the  parameter  is  made  by  decreasing  the  torque 
constant  ratio  by  20%.  The  nonlinear  finite  element  model  is  integrated  to  the 
measurement  times  and  the  resulting  model  values  are  used  in  Eq.  (160)  to  update  r.  Table 
22  sxunmarizes  the  results  of  the  recovered  torque  constant  ratio! 


Table  22.  Recovery  of  torque  constant  ratio 

True  value 

1.0000 

Starting  Estimate 

0.8000 

Value  after  1  iteration 

0.9872 

Converged  value 

0.9999 

The  ratio  recovered  after  one  iteration  is  very  close  to  the  true  value  and  is  essentially  the 
same  after  two  iterations.  The  ‘true’  torque  constant  is  recovered  quickly  and  accurately 
when  measurement  errors  are  neglected. 

Measurement  Errors 

Errors  are  introduced  into  the  simulated  measurements  using  a  Gaussian  random 
number  generator  and  an  associated  standard  deviation  for  each  value.  These  standard 
deviations  range  from  0.0008  to  0.0020  for  the  w,  and  m,  measurements.  The  parameter 
recovery  program  is  repeated  using  the  weighted  measurements.  Table  23  shows  the 
resulting  torque  constant  ratio  recovered. 


1  Table  23.  Recovery  of  ratio  with  measurement  errors 

True  values 

1.0000 

Initial  estimate 

0.8000 

Recovered  value 

0.9567 

Recovered  a 

0.0514 

The  torque  constant  ratio  is  recovered  after  two  iterations  and  is  within  4%  of  the  true 
value.  The  standard  deviation  for  this  parameter  is  small  and  indicates  this  ratio  is  well 
determined. 


CHAPTER  VI 


CONCLUSIONS 

In  this  dissertation,  several  nonlinear  static  and  dynamical  models  for  large 
deformations  of  flexible  body  chains  have  been  investigated.  A  novel  formulation  is 
introduced  for  nonlinear  beam  static  deformations  using  arc  length  as  the  independent 
spatial  variable.  This  formulation  is  validated  by  comparison  to  several  published 
solutions  and  experiments.  For  transient  dynamic  response,  a  geometrically  exact  finite 
element  model  is  developed  fi’om  a  static  model  proposed  by  Epstein  and  Murray. 

The  focus  of  this  dissertation  is  the  development  of  a  three  step  approach  to  system 
identification: 

1 .  A  series  of  static  deformation  experiments  are  used  to  recover  the  parameters 
affecting  the  potential  energy  stored  by  the  beams. 

2.  The  beam  system  is  released  from  the  statically  deformed  position  and  free 
response  measurements  are  used  to  recover  the  inertia  type  parameters. 

3.  A  forced  response  experiment  is  used  to  recover  actuator  model  parameters. 
Measurement  and  a  priori  covariance  estimates  can  be  propagated  through  the 
sequence  of  parameter  recoveries. 

Based  on  the  results  of  tins  approach,  the  following  conclusions  are  drawn. 

1 .  The  beam  parameters  affecting  the  potential  energy  of  the  system  can  be  quickly 
and  accurately  recovered  using  the  arc  length  model  and  a  family  of 
deformation  shapes  from  two  types  of  static  experiments.  The  stiffness 
coefficients  are  well  determined  in  the  presence  of  measurement  errors  for 
appropriately  designed  experiments.  The  stiffiiess  coefficient  of  the  first  beam 
is  better  determined  since  the  portion  of  the  system  closest  to  the  clamped 
boimdary  condition  deforms  the  most  and  parameters  affecting  that  response  are 
more  easily  estimated. 


2.  The  beam  parameters  affecting  the  free  response  of  the  system  can  be  recovered 
using  a  nonlinear  finite  element  model  and  the  observed  motion  due  to  the 
release  of  the  constraining  wire  from  the  static  experiment.  The  coefficients  of 
the  friction  forces  between  the  tip  and  elbow  masses  and  the  supporting  table 
are  well  determined  with  and  without  the  presence  of  measurement  errors.  The 
mass  density  per  unit  length  for  each  beam  can  be  determined  assuming  perfect 
measurements  but  are  difficult  to  recover  in  the  presence  of  measurement  errors 
unless  large  external  (e.g.  friction)  forces  are  present.  The  mass  densities  do  not 
have  a  large  affect  on  the  free,  lightly  damped  motion  of  the  system  and  are 
difficult  parameters  to  determine.  The  parameters  in  the  second  beam  are  better 
estimated  than  the  parameters  in  the  first  beam  for  this  clamped/free  system  due 
to  larger  dynamic  motion  at  the  tip  of  the  system. 

3.  The  motor  parameters  affecting  the  forced  response  of  the  system  can  be 
recovered  using  the  nonlinear  finite  element  model  and  the  measured  response 
over  a  short  transient.  For  the  cases  studied,  motor  torque  constants  are  well 
determined  in  the  presence  of  realistic  measurement  errors. 

4.  The  sequential  process  permits  the  dimensionality  of  the  estimation  process  to 
be  increased  as  more  information  is  available  thus  enhancing  convergence  in 
high-dimensional  parameter  estimation  processes. 

Several  difficulties  were  encountered  and  I  suggest  the  following  issues  for  further 
investigation: 

1 .  The  arc  length  approach  works  well  for  the  nonlinear  static  analysis  but  cannot 
be  easily  extended  to  the  dynamic  analysis  under  general  modeling 
assumptions.  The  nonlinear  finite  element  model  could  be  used  to  recover  the 
potential  energy  parameters  for  model  consistency  in  bringing  forward  die 
recovered  parameters  and  their  associated  uncertainties. 

2.  The  nonlinear  finite  element  model  works  well  for  simulating  the  dynamic 
response  and  recovering  the  free  and  forced  vibration  parameters.  However,  a 
small  time  step  must  be  used  for  the  stability  of  the  numerical  integration  which 


leads  to  a  large  amount  of  computer  time  for  tracking  the  motion  over  a  long 
time  period.  Optimization  of  the  computer  code  and  variable  time  steps  could 
be  used  to  shorten  the  amount  of  computer  time  needed.  Some  new  method  is 
needed  to  do  order  reduction  for  this  class  of  problems. 

3.  A  method  for  carrying  forward  the  parameter  estimates  and  uncertainties  in  the 
sequence  of  parameter  recovery  experiments  has  been  developed.  Simulation 
studies  of  this  process  at  each  step  needs  to  be  further  investigated. 
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